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ABSTRACT 

This  thesis  is  a  continuation  of  previous  work  to  apply 
modern  multiparameter  estimation  techniques  to  the  problem 
of  estimating  attrition  rates  for  a  large  number  of  small 
inventory  cells  in  manpower  planning  models  used  by  the  U.S. 
Marine  Corps.  The  main  advances  involve  the  promising 
introduction  of  empirical  Bayes  (non-constant  shrinkage) 
techniques,  recognition  of  the  non  symmetric  nature  of  the 
errors  with  a  response  to  this,  and  some  insight  into  all 
aggregation  plans  that  should  help  provide  greater  stability 
for  the  estimation  methods.  In  addition,  the  roles  of  some 
middle  level  methodological  choices  are  explored. 
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THESIS  DISCLAIMER 


The  reader  is  cautioned  that  computer  programs 
developed  in  this  research  may  not  have  been  exercised  for 
all  cases  of  interest.  While  every  effort  has  been  made, 
within  the  time  available,  to  ensure  that  the  programs  are 
free  of  computational  and  logic  errors,  they  cannot  be 
considered  validated.  Any  application  of  these  programs 
without  additional  verification  is  at  the  risk  of  the  user. 
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I.   INTRODUCTION 

A.   PURPOSE 

In  support  of  project  OPUS  (Officer  Planning  and  Utility 
System) ,  Headquarters,  U.S.  Marine  Corps,  in  1985,  requested 
assistance  in  exploring  new  methods  of  generating  manpower 
loss  rates  in  order  to  improve  upon  the  currently  used 
method  of  force  rate  analysis.  In  response,  students  and 
faculty  in  the  Operations  Research  department  at  the  Naval 
Postgraduate  School  began  exploring  the  use  of  modern 
multiparameter  estimation  techniques  for  this  problem. 
Special  emphasis  is  placed  on  the  small  cell  problem,  i.e., 
categories  of  officer  skill,  grade  and  length  of  service 
which  have  low  inventory  figures.  Historically,  rate 
estimators  for  small  cells  are  unstable  and  a  large  number 
of  these  cells  exist. 

This  paper  builds  upon  previous  work  on  this  problem. 
The  main  advances  are  to  append  new  measures  of 
effectiveness  as  requested  by  the  Navy  Personnel  Research 
and  Development  Center  (NPRDC),  introduce  a  class  of 
empirical  Bayes  estimators,  and  explore  the  effects  of  some 
middle  level  choices  in  applying  the  new  and  existing 
methods.  The  ultimate  goal  is  to  refine  the  techniques 
presented  here  in  order  to  validate  a  clear  policy  for 
predicting  loss  rates. 


B .   BACKGROUND 

For  an  introduction  into  Marine  Corps  policy  concerning 
manpower  planning,  the  reader  is  referred  to  a  thesis 
submitted  in  September  1985,  at  the  Naval  Postgraduate 
School,  by  Major  D.D.  Tucker  [Ref.  1],  who  presents  a 
detailed  background  into  the  Marine  Corps  officer  structure 
and  the  manpower  planning  process.  Tucker  also  provides 
basic  attrition  rate  theory  and  calculates  overall  attrition 
rates  in  several  different  formats  as  they  pertain  to  the 
Marine  Corps.  As  the  aggregation  rate  begins  to  grow 
smaller  with  refinement,  Tucker  illustrates  the  irregular 
behavior  of  losses  due  to  voluntary  attrition  and  the  small 
numbers  of  losses  obtained  in  such  aggregates.  Tucker  [Ref. 
l:p.  50]  introduces  the  James-Stein  technique  of  loss  rate 
estimation  as  a  shrinking  of  individual  cell  averages  toward 
a  grand  average  in  order  to  reduce  the  risk,  or  squared 
differences  of  forecast  and  actual  values,  the  goal  being  an 
improvement  over  the  classical  maximum  likelihood  estimator, 
referred  to  from  now  on  as  the  MLE .  Tucker  chose  the  ranks 
of  First  Lieutenant  and  Lieutenant  Colonel  for  evaluation 
using  military  occupational  specialty  (MOS)  groups  of  combat 
support  and  ground  combat  (three  MOS ' s  were  selected  for 
each  group) .  These  aggregation  schemes  were  carried  through 
other  studies  and  are  included  for  historical  purposes  in 
this  paper. 


Another  thesis  submitted  at  the  Naval  Postgraduate 
School  in  March  1986,  by  Major  J.R.  Robinson  [Ref .  2] 
focused  on  a  technique  called  limited  translation  of  the 
James-Stein  estimator.  This  particular  model  attempts  to 
minimize  the  risk  mentioned  above  by  reducing  the  shrinkage 
of  rate  estimates  for  cells  which  fall  outside  a  certain 
range  of  values  centered  on  the  grand  mean.  It  was  thought 
that  shrinking  all  cells  by  the  same  rate  toward  this  grand 
mean  may  be  unwise  since  the  attrition  rates  for  those 
particular  cells,  those  farthest  from  the  grand  mean,  may 
contain  attributes  which  are  different  than  the  majority  of 
the  cells  in  the  aggregate.  Additionally,  Robinson's  work 
carries  forward  the  study  done  by  Tucker  concerning  MLE  and 
James-Stein  estimators. 

These  feasibility  studies  were  limited  in  scope  due  to 
the  format  of  the  summary  date  tape  available  (acquired  from 
NPRDC) ,  which  allows  only  coarse  cell  definition.  A  cell  is 
defined  to  be  a  cross-classification  of  the  forty  military 
occupational  fields  (OF) ,  thirty-one  lengths  of  service 
(LOS) ,  and  ten  grades  for  a  total  of  12,400  categories  for 
manpower  planning  purposes  [Ref.  2:p.  10].  The  data  does 
not  distinguish  between  limited  duty  officers,  those 
officers  specifically  designated  for  limited  duty  within 
certain  MOS ' s ,  and  unrestricted  officers  (regular  and 
reserve)  [Ref.  l:p.  20].  These  factors  dictated  a  broad 
aggregation  scheme   that  created  stable  cell  inventories  and 


allowed  the  common  variance  assumption  crucial  to  their 
work.  But  some  of  the  application  models  require  more 
refined  information  along  those  lines  and  experiments  with 
finer  aggregation  levels,  necessary  in  some  real  cases, 
caused  this  assumption  to  be  violated. 

Tucker  and  Robinson  did  conclude  that  the  current  scheme 
of  attrition  rate  prediction  could  be  improved  upon  by  their 
methods,  although  no  dominant  scheme  was  uncovered.  Tucker 
and  Robinson  did  identify  problem  areas  in  their  studies 
that  required  additional  work.  Tucker  [Ref .  l:p.  71]  cited 
a  need  for  a  better  aggregation  method  to  produce  a  more 
homogeneous  attrition  behavior  in  all  cells.  Also,  small 
probabilities  of  loss  within  a  cell  were  not  dealt  with 
successfully.  Robinson  [Ref.  2:p.  32]  revealed  the  small 
cell  problem  and  the  inability  to  normalize  the  cell  means 
or  stabilize  the  variance  with  a  data  transformation.  The 
recommendations  of  the  previous  studies  highlighted  the  need 
for  a  more  refined  data  tape  with  current  information  to 
include,  among  others,  full  Military  Occupational  Specialty 
(MOS)  information,  grade  separation  to  include 
regular/reserve  status,  promotion  zone  data  and  breakdown  of 
attritions  by  type. 

C.   PROGRESS 

The  project  of  loss  rate  estimation,  sponsored  by  MPRDC , 
is  currently  moving  into  the  implementation  stage.  An 
operational   data   tape   is   now   available   which   includes 
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detailed  information  on  Marine  Corps  officers  for  the  years 
1977-86.  Work  has  begun  at  the  Naval  Postgraduate  School  to 
break  out  the  data  in  a  useful  format  and  its  availability 
in  now  nearing  completion. 

In  addition  to  the  results  stated  earlier,  Tucker  and 
Robinson  provided  loss  rate  estimates  for  small  cells  with 
no  attrition,  i.e.,  MLE  equal  to  zero,  and  risk  effects  of 
the  different  schemes  studied.  These  results  provided  the 
springboard  for  this  study  and  the  requirements  for 
additional  future  studies. 

Captain  R.W.  Larsen  [Ref.  3]  has  developed  a  promising 
aggregation  scheme  based  on  cluster  analysis,  i.e.,  a 
classification  scheme  to  aggregate  cells  which  reflect  a 
greater  degree  of  homogeneity  in  attrition  rates  than  are 
allowed  using  current  aggregation  methods.  Table  1  displays 
the  resulting  scheme,  given  in  year  of  current  service  (YCS) 
groups  proposed  for  implementation  by  Larsen. 


Table  1.     AGGREGATION  METHOD  PROPOSED  BY  LARSEN 


MPS  Category  Bounded  YCS  Groups 

Fixed-Wing  Pilots  (1-6,8-19)    (7)      (20-25)  (26) 

Rotary- Wing  Pilots  (1-5,8-19)    (6,7)    (20-25)  (26) 

Naval  Flight  Officers  (1-5,8-19)    (6,7)   (20-25)  (26) 

Lawyers  (1-6,8-19)    (7)      (20-25)  (26) 

All  Else  -  (1-3,6-19)    (4,5)   (20-25)  (26) 


The  parentheses  encompass  YCS  groups  which  behave  similarly. 
In  each  case  the  grade  is  fixed,  i.e.,  one  should  aggregate 
over  YCS  before  aggregating  over  grade. 

This  information  will  allow  for  more  successful 
application  of  the  loss  rate  estimation  techniques  of 
interest  here  when  the  new  data  tape  is  provided. 

B.   GOALS 

In  an  attempt  to  refine  the  methods  of  attrition  rate 
estimation  presented  in  the  previous  pilot  stud!^b,  the 
following  goals  were  set  for  this  paper: 

1.  Modify  the  existing  methods  to: 

a)  use  the  refined  data  format, 

b)  extend  the   work  to   the  case  of  unequal  or  non 
constant  variance  in  the  small  cell  inventory 
problem, 

c)  study   alternate   transformation   inversion 
techniques , 

d)  introduce  and  evaluate  additional   measures  of 
effectiveness  suggested   by  NPRDC  to  include,  but 
not  limited  to,  cross  validation  underages   and 
overages  of   the  deviation  between  forecasts  and 
actuals , 

e)  examine  the  graphical  effect  of  shrinkage   in  the 
original  scale. 

2.  Introduce  empirical  Bayes  estimation  method  in  several 
forms  for  attrition  rates  for  consideration  as  another 
option  to  solve  the  problem. 


II.   EXTENSIONS  TO  ESTIMATION  METHODS 

A.  GENERAL 

As  stated  in  Chapter  I,  Tucker  and  Robinson  published 
pilot  studies  introducing  the  James-Stein,  limited 
translation  James-Stein,  Maximum  Likelihood  and  Transform 
Scale  Cell  Average  estimators  as  alternatives  to  the 
aggregate  methods  previously  used  for  manpower  planning  in 
the  Marine  Corps.  Their  performance  was  promising,  and  this 
chapter  will  investigate  ways  to  enhance  their  performance 
and  possibly  develop  a  dominant  estimator  for  future 
attrition  prediction.  Additionally,  NPRDC  has  expressed 
interest  in  additional  measures  of  effectiveness  for  these 
estimators  which  may  lead  to  a  sharper  direction  for  the 
goal  of  producing  a  viable  estimation  policy.  This  chapter 
concentrates  on  the  following  issues  concerning  previous 
studies . 

B.  NON  CONSTANT  VARIANCE 

For  a  cell  having  an  inventory  n,  and   attrition  rate  p, 
the   number   of   attritions   y   is  described  by  the  Binomial 
(n,p)  distribution.   The  variance  of  the  estimator  of 
p  =  y/n  is  given  by  the  familiar  formula 

var(p)  =Pd  ~  P> 

y  n  (2.1) 


The  Freeman-Tukey  double  arcsine  transformation  was  used 
by  Tucker  [Ref.  l:p.  55]  and  Robinson  [Ref.  2:pp.  74-75]  in 
order  to  map  the  raw  losses  to  the  transformed  scale.  This 
formula  is  given  by 

x  =  0.5[n  +  0.5]  ^  2  (sin-1  [, — ^H 1]         '2-2^ 

(n  +  1 ) 

Note  that  n=n  ( t )  may  change  with  time  and  thus,  so  would  x. 
The  question  of  averaging  over  time  before  or  after  applying 
the  Freeman-Tukey  transformation  is  not  discussed  in  this 
paper.  One  of  the  goals  of  using  this  transformation  is  to 
stabilize  the  variance  of  this  estimate.  In  the  work  of 
Tucker  and  Robinson,  the  assumption  of  normally  distributed 
random  variables  with  common  variance  was  the  setting  in 
which  the  James-Stein  estimator  was  expected  to  perform 
well.  Indeed,  for  n  moderate  or  large  and  p  not  too 
extreme,  the  variance  of  the  Freeman-Tukey  transformation  is 
approximately  one.  But  as  discovered  by  Robinson,  this 
assumption  does  not  hold  well  when  applied  to  cells  with 
small  inventory.  A  more  careful  look  at  the  variance  of  the 
transformed  data  is  taken  here,  with  particular  interest  in 
the  region  of  unstable  variance  where  n  and  p  are  small. 

Some  exploratory  graphical  work  has  revealed  that  the 
natural  log  of  the  variance  of  x,  the  transformed  value, 


behaves  quite  well  as  a  linear  function  of  u(  =E[X])  and  u- 
1  in  this  unstable  zone.  This  leads  to  the  interpolatory 
formula : 

var(x)  =  1.683511- ■  8934  (u  -  l)-9981,  (2.3) 

1.001  <;  u  ^  2.2 

The  upper  limit,  u  =  2.2,  is  the  value  for  which  var(x)  =  1. 
Since  the  least  squares  fitting  process  was  applied  to  the 
logarithms,  u  can  not  be  allowed  to  fall  to  one  or  below. 
The  lower  limit,  u  =  1.001,  is  an  arbitrary  value  that  meets 
this  constraint.  The  coefficients  were  validated  by  a 
linear  regression  on  the  natural  log  of  this  formula.  The 
actual  variance  function  for  n  =  1,2,3,4,5,7,10  appear  in 
Figure  1,  and  the  fit  is  remarkably  good  for  n  >  3.  The 
upper  graph  shows  the  region  of  special  interest  for  this 
paper  and  this  portion  of  the  curve  is  fitted  by  the  formula 
in  Equation  2.3  above.  The  lower  plot  exhibits  the  outer 
tale  effect  yet  to  be  explored.  The  small  values  of  n(<3) 
do  not  maintain  the  stable  variance  region,  i.e.,  var[x]=l, 
for  long  before  falling  back.  As  the  value  for  n  increases, 
the  more  stable  the  variance  becomes.  Table  2  is  a  partial 
display  of  the  residuals  for  selected  n  values,  fifteen 
equally  spaced  values  in  the  range  1.001  <  u  ^  2.2,  computed 
by  the  linear  regression  formula  above.  For  values  of  n  >  2 
the  difference  between  actual  and  fitted  values  are 
considered  acceptable  for  this  work. 
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Table  2.     RESIDUAL  LISTING  FOR  SELECTED  VALUES  OF  N 

N-l 

N-2 

N  =  4 

N-5 

N=10 

■ 

0.07191 

0.03912 

0.02139 

0.017 

0.008313 

0.08466 

0.02569 

0.002853 

-0.001695 

-0.01106 

0.1033 

0.0239 

-0.005448 

-0.01049 

-0.0206 

0.122 

0.03002 

-0.00655 

-0.01 1S5 

-0.02257 

0.1385 

0.04068 

-0.004026 

-0.009566 

-0.02067 

0.15 

0.053S8 

0.0000676 

-0.005727 

-0.0166 

0.1552 

0.06853 

0.004027 

-0.002103 

-0.01291 

0.1527 

0.0S234 

0.007447 

0.0009541 

-0.009981 

0.1414 

0.09481 

0.00948 

0.002557 

-0.008299 

0.1209 

0.1055 

0.009825 

0.002448 

-0.008128 

0.08958 

0.1129 

0.008369 

0.0004583 

-0.009929 

0.04832 

0.1169 

0.004997 

-0.003511 

-0.01327 

-0.005254 

0.1171 

0.0001116 

-O.OOS93 

-0.0182 

-0.06914 

0.113 

-0.006798 

-0.01654 

-0.02463 

-0.1466 

0.104 

-0.01534 

-0.025S5 

-0.03366 

The  boundaries  for  u,  in  the  unstable  region,  as  a 
function  of  n  and  p  are  provided  in  Figure  2.  The  formula 
derived  was  quite  adequate  for  this  study. 


C.   ALTERNATE  TRANSFORM  INVERSION  METHODS 

There  are  three  issues  in  the  inversion  of  the  data  from 
the  transformed  space.  The  first  is  that  of  selecting 
values  of  n,  i.e.,  an  average  of  the  n(t)'s  appearing  in  the 
estimation  year  inventories.  This  issue  is  discussed  later 
in  the  chapter,  with  the  development  of  a  simulation  for 
comparing  mean  inventory  values.  The  second  is  that  of  the 
inversion  of  the  sum  of  the  two  arcsine  functions.  The 
third  involves  the  case  of  time  averaging  in  the  transformed 
scale  and  the  question  of  what  matching  form  of  average 
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inventory  should  be  used  in  the  inversion.  The  third  issue 
is  not  covered  in  this  paper,  and  is  still  pending. 
Treatment  of  the  second  issue  will  be  provided  first. 

The  method  used  in  the  pilot  studies  of  Tucker  and 
Robinson,  and  the  most  transparent  transform  inversion 
formula  is: 

Pl  =  0.5[1  -  sin((n  -    1/2)1/2)]  (2.4) 

This  will  be  referred  to  as  the  Basic  Inversion.  The 
formula  offered  by  Carter  and  Rolph  [Ref.  4]  is  given  by: 


a   -  ,--;~2  t^    \    >  (1  ~  2sin2  9i  )  ,2 
pt  -  smMGi)  +  I(4(n/k)  +  2)   IBi 


2.5) 


where  9i  is  x/ (n  +  1/2) 1/2  and  Bt  is  the  estimated  amount  of 
shrinkage  in  the  transformed  scale  depending  on  the  Bayes 
model  used  (constant  or  proportional  prior) ,  k  is  the  number 
of  cells  and  n  is  inventory. 

Apparently,  Carter   and  Rolph   invert  prior  to  shrinking 
and  the  above  formula  provides  a   shrinkage  in   the  original 
scale.   Rewriting  the  above,  it  appears  to  shrink  towards 
f>  =  1/2,  and  is  given  by: 

»'  "  1/2(4!n/k)'+  2>  +  ■*»»'<*  XI  ~« (n/kf  +  21       (2-6) 

An  exact  Freeman-Tukey  inversion  formula  developed  by 
Miller  [Ref.  5]  was  recently  uncovered.  It  is  exact  in  the 
sense  that  for  a  single  cell,  the  formula  given  by: 
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p(t)  =  0 . 5  f 1  -  sgn(cos  t)  [1  -  (sin  t  + 
sin  5  -  1/sin  t 


2 


n 


I1 /2  !  (2.7) 


where  t  =  x/(n  +  1/2)1/2  +  n/2  and  a  is  the  average  central 
inventory,  returns  the  empirical  rate  y/n.  This  will  be 
referred  to  as  the  Freeman-Tukey  Exact  (FTE)  inversion. 
Since  this  formula  is  applied  using  an  average  value  of  n, 
there  may  be  concern  that  it  oscillates  between  successive 
integer  values  of  n.  Figure  3  indicates,  using  values  of 
n=5  and  10,  that  it  is  quite  smooth  as  an  interpolation 
formula.  In  the  article,  Miller  suggested  using  the 
harmonic  mean  of  the  n  values  provided  by  the  time  changes . 
This  issue  is  examined  later  in  this  chapter.  The  Miller 
formula  was  incorporated  in  the  Marine  Corps  Inventory 
Attrition  Analysis  Program  provided  in  Appendix  B. 

D.   SHRINKAGE  VISUALIZATION 

When  shrinkage  is  a  fixed  amount  for  all  cells,  it  is 
readily  visualized,  for  all  degrees  of  shrinkage,  by  the 
linear  diagrams  presented  in  the  works  of  Tucker  [Ref.  l:p. 
52]  and  Casella  [Ref.  6].  Robinson  [Ref.  2:p.  20]  also 
presented  diagrams  of  the  limited  translation  option  and  how 
the  shrinkage  is  affected.  All  of  these  diagrams  refer  to 
the  transformed  scale. 

Of  more  immediate  interest  is  the  question  of  how  to 
visualize  shrinkage  diagrams  of  this  type  in  the  original 
scale.   Illustrations  of  diagrams  of  this  type  are  provided 
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using  both  the  basic  inversion  formula  used  by  both  Tucker 
and  Robinson,  and  the  FTE  formula.  Figure  4  shows  how  the 
diagram  by  Casella  [Ref.  6:p.  84]  looks  in  the  original 
scale  using  the  same  value  of  n  for  all  cells.  An 
interesting  note  to  consider  is  the  effect  of  the  difference 
between  inversion  formulas,  i.e.,  in  the  graph  which 
combines  the  Basic  and  FTE  inverses,  the  FTE  inversion  tends 
to  shrink  at  a  slower  rate,  as  values  get  further  from  the 
mean.  For  a  modest  transition  into  the  case  of  non  constant 
n,  Figure  5  provides  for  two  values  of  n.  Finally,  Figure  5 
illustrates  the  point  with  many  different  values  of  n,  taken 
from  actual  Marine  Corps  officer  inventory  data.  In  this 
case,  note  that  Basic  inversion  shrinks  to  one  original 
scale  value  while  the  FTE  inverse  is  not  as  focused. 

E.   CHOICE  OF  AVERAGE  INVENTORY  VALUES  OVER  TIME 

A  Fortran  simulation  was  developed,  included  as  Appendix 
C,  to  test  the  performance  of  each  of  several  common  means 
(arithmetic,  geometric  and  harmonic)  in  the  comparison  of 
the  inverted  estimates  with  the  empirical  y/n.  Values  of  n 
were  produced  from  a  Poisson  distribution  with  rate  X  (1  <  ,\ 
<  20) ,  and  attrition  values,  y,  were  produced  from  the 
Binomial  (n,p)  distribution  where  a  was  in  the  range  .01  S  p 
>.40.  The  n  and  y  values  were  then  transformed  by  using  the 
Freeman-Tukey  formula.  An  exhaustive  study  of  inversion 
variations  with  both  the  FTE  formula  and  the  basic  formula 
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used  by  Robinson  were  compared.  Table  3  contains  selected 
outputs  which  are  typical  of  the  broader  study  and  shows  how 
the  means  performed  against  each  other  in  three  different 
computations  listed  below. 

1.  The  BASIC  Method  utilizes  the  Basic  inversion  formula 
and  the  average  transformed  scale  value. 

2.  The  FTE1  Method  utilizes  Equation  2.7  with  the  t 
values  equal  to  the  average  (over  replications,  i.e., 
time)  of  x/(n  +  1/2) 1/2  and  >:  given  by  Equation  2.2. 
More  explicitly,  see  Equation  A. 16. 

3.  The  FTE2  Method  utilize  Equation  2.7  with  the  t 
values  equal  to  the  average  (over  replication,  i.e.. 
time)  of  x  directly,  and  x  given  by  Equation  2.2 
followed  with  a  division  by  (n  +  1/2)1/2 

In  the  context  of  this  study,  the  arithmetic  mean  tended 

to  perform  as  well  or  better  than  the   geometric  or  harmonic 

means.    The  choice   of  arithmetic  mean  also  circumvents  the 

problem  of  an  undefined  term   when   nt   =   0   is  encountered 

while  using  the  harmonic  of  geometric  means. 

F.   QUESTIONS  RELATING  TO  VALIDATION 

The  pilot  studies  by  Tucker  and  Robinson  both  used  four 
years  (1977-80)  for  estimation  of  the  attrition  rates  and 
the  remaining  three  years  (1981-83)  for  validation.  The 
estimation  rates  were  used  in  an  attempt  to  forecast  the 
values  for  the  three  validation  years .  The  quality  of  the 
forecasts  deteriorates  as  lead  time  increases  in  the  future. 
This  was  anticipated,  in  fact,  Rowe  et  al . ,  at  NPRDC  have 
shown  a  time  series  effect  in  their  work  forecasting 
attritions  in  the  Navy  [Ref.  7]  and  are  currently 
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Table 

3.     SIMULATION    OUTPUT    OF   MEAN 

VALUES    FOR    INVERSION 

UTILIZATION 

X_ 

REPS 

P 

METHOD 

ARITHMETIC 

GEOMETRIC    HARMONIC 

2.0 

40 

0.05 

BASIC 

0.1068 

0.1176 

0.1308 

FTE1 

0.0252 

0.01S9 

0.0122 

FTE2 

0.0181 

0.0208 

0.0245 

2.0 

40 

0.10 

BASIC 

0.1399 

0.1537 

0.1708 

FTE1 

0.0535 

0.0458, 

0.0369 

FTE2 

0.0489 

0.0558 

0.0648 

2.0 

40 

0.20 

BASIC 

0.2402 

0.2629 

0.2906 

FTE1 

0.1900 

0.1829 

0.1739 

FTE2 

0.1657 

0.1869 

0.2144 

2.0 

40 

0.40 

BASIC 

0.3713 

0.4038 

0.4429 

FTE1 

0.3442 

0.3403 

0.3353 

FTE2 

0.3325 

0.3715 

0.4213 

4.0 

40 

0.05 

BASIC 

0.0750 

0.0815 

0.0921 

FTE1 

0.0235 

0.0200 

0.0147 

FTE2 

0.0200 

0.0221 

0.0258 

4.0 

40 

0.10 

BASIC 

0.1280 

0.1387 

0.1564 

FTE1 

0.0768 

0.0724 

0.0655 

FTE2 

0.0726 

0.0798 

0.0922 

4.0 

40 

0.20 

BASIC 

0.1965 

0.2125 

0.2386 

FTE1 

0.1446 

0.1406 

0.1339 

FTE2 

0.1485 

0.1628 

0.1S69 

4.0 

40 

0.40 

BASIC 

0.4382 

0.4694 

0.5185 

FTE1 

0.4564 

0.4558 

0.4549 

FTE2 

0.42S0 

0.4638 

0.5223 

8.0 

40 

0.05 

BASIC 

0.0552 

0.05S0 

0.0616 

FTE1 

0.0181 

0.0166 

0.0148 

FTE2 

0.0156 

0.0165 

0.0177 

8.0 

40 

0.10 

BASIC 

0.1226 

0.1285 

0.1363 

FTE1 

0.0922 

0.0904 

o.osso 

FTE2 

0.0828 

0.0873 

0.0933 

8.0 

40 

0.20 

BASIC 

0.2035 

0.2131 

0.2255 

FTE1 

0.1778 

0.1762 

0.1741 

FTE2 

0.1708 

0.1797 

0.1915 

8.0 

40 

0.40 

BASIC 

0.3799 

0.3962 

0.4170 

FTE1 

0.3700 

0.3693 

0.3684 

FTE2 

0.3662 

0.3837 

0.4065 
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recommending  a  second  order  auto  regression  scheme  for  this 
purpose.  In  this  regard,  the  work  of  Tucker  and  Robinson 
with  validation  may  be  viewed  as  forecasting  by  persistence , 
and  in  light  of  the  NPRDC  work,  can  hardly  be  expected  to 
perform  well  in  an  absolute  sense.  Their  work  does, 
however,  provide  some  comparisons  that  can  be  useful  in 
choosing  among  alternative  techniques,  and  much  can  be 
learned  from  them. 

The  first  measure  considered  is  the  Mean  Squared  Error 
(MSE)  of  forecasts.  The  shrinkage  estimation  theory 
provides  that  these  values  should  be  about  unity  if 
persistence,  independence  over  time,  was  agreeable.  That 
is,  useful  validations  can  be  obtained  without  modeling  the 
time  dependence  of  the  attrition  process.,  Only  for  a  very 
few  cases  does  this  hold  (i.e.,  Warrant  Officers,  see 
Tucker).  However,  if  one  method  of  estimation  consistently 
produces  noticeably  smaller  MSE  values  chan  do  the  other 
methods,  this  would  provide  sufficient  support  for  that 
method's  continued  development. 

The  second  measure  is  designed  to  consider  performance 
in  the  original  scale  and  is  patterned  after  a  chi-square 
statistic 

-,   (actuals  -  forecasts)2  ,0  ox 


j        forecasts 

If  the  k  cells  are   independent,   then   equation   2.3  is 
approximately  a  chi-square  random  variable  with  k  degrees  of 
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freedom.  Hence,  the  man  value  is  about  k  and  the  standard 
deviation  is  (2k)1/2.  Once  again,  we  generally  do  not  see 
values  this  small,  but  use  it  for  comparison  purposes.  As  a 
side  note,  the  value  of  k  when  using  the  MLE  is  generally 
smaller  than  that  of  the  other  estimators  due  to  the 
omission  of  those  cells  with  MLE  equal  to  zero  over  all 
estimation  years,  i.e.,  zero  forecasts.  This  provides  and 
unnatural  advantage  for  the  MLE,  which  must  be  kept  in  mind 
when  interpreting  the  results. 

Two  disturbing  features  have  emerged  from  the 
exploratory  analysis.  Often  a  technique  that  performs  well 
for  some  cases  using  the  MSE  measure  is  not  comparably 
supported  by  the  chi-square  measure.  It  seems  the  two 
measures  do  not  weight  the  common  features  equally,  and 
further  study  is  necessary  to  develop  understanding  of  the 
results. 

Discussion  with  personnel  at  NPRDC  concerning  these 
matters  has  led  to  the  introduction  of  a  third  measure, 
which  will  be  called  the  Mean  Absolute  Deviation  (MAD). 
This  value  is  computed  in  the  original  scale.  Moreover, 
since  the  cost  structure  of  over-forecasting  is  not 
symmetric,  we  separate  the  two  parts  into  the  average 
overage   deviation  given  by: 


i  k 

ODEV  =  ±  S  [forecasts  -  actuals]*       (2.9) 
k   1 
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and  average  underage   deviation 


i   k 

UDEV  =  r-  2Z  [actuals  -  forecasts]*      (2.10 
K      1 


where   the   "+"   refers   to   the   positive  part,  as  computed 
separately.   Finally, 

MAD  =  ODEV  +  UDEV  (2.11) 

This  measure  also  has  shown  that  in  most  cases,  the 
overage  exceeds  the  underage  value,  which  in  turn  may  lead 
to  forecasting  the  need  for  an  excess  number  of  Marine 
officers. 
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III.   EMPIRICAL  BAYES  ESTIMATION 

A.  BACKGROUND 

Although  the  studies  of  Tucker  and  Robinson  showed 
improvement  over  the  current  methods,  no  clear  procedure  was 
established  and  no  plans  for  implementation  were  drafted. 
Problem  areas  were  identified  while  working  with  small  cell 
inventories  and  low  loss  rates.  In  the  transformed  scale, 
the  James-Stein  assumption  of  normal  distribution  of  the 
observations  with  equal  variance  did  not  hold  well  in  such 
cases.  An  alternate  method  of  loss  rate  estimation  is 
required  to  deal  with  this  special  case. 

B .  INTRODUCTION 

Empirical  Bayes  is  simply  a  generalization  of  the  James- 
Stein  estimator  for  use  with  small  cell  inventories  and 
unstable  variance.  Random  variables  are  assumed  to  be  of 
the  form: 

Xi  ~  N(8i  ,di  )  ,  i=l,2 k  (3.1) 

where  9i  are  the  unknown  parameters  to  be  estimated  and  di 
are  known  variances.  The  Bayesian  assumption  is  that  the 
0i  's  are  themselves  from  some  distribution,  in  this  case 
also  from  the  normal  distribution  with  mean,  v,  and 
variance,  pi  a.  The  proportionality  constant  pi  is 
normalized  so  that 
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S  /Oi  =  k,  i=lP2 k  (3.2) 

where  k  is  the  number  of  observations  or  cells.   Thus, 

9i  "  N(v,^ia),  i=l,2 k.  (3.3) 

This  distribution  is  known  as  the  prior  distribution  of  Oi 
and  represents  the  additional  assumed  information  on  the 
value  of  9  before  observing  X  whose  distribution  depends  on 
0  [Ref.  8:p.  553].  The  use  of  this  information  produces  an 
estimate  which  is  a  weighted  average  of  the  prior  mean,  v, 
and  the  sample  estimate  x.  In  the  empirical  Bayes  context, 
the  values  of  v  and  a,  are  also  estimated  from  the  data 
since  all  the  information  necessary  is  contained  in  the 
marginal  distribution  of  Xi  (unconditional  on  ©i )  [Ref.  6:p. 
83] .   This  chapter  focuses  on  that  method. 

For  a  more  detailed  explanation  of  the  Bayesian  family 
of  estimators,  the  reader  is  referred  to  Casella  [Ref.  6] , 
Carter  and  Rolph  [Ref.  4]  and  Efron  and  Morris  [Ref.  9] . 

C.   THE  ESTIMATION  METHOD 

Fay  and  Herriot  [Ref.  10]  discuss  the  empirical  Bayes 
technique  and  provide  an  estimate  for  the  prior  mean  as 


S   xi  /(a*  +  di  ) 

vk    =  _J ,  i=l,2 k  (3.4) 

ZZ   l/(a*  +  di  ) 

i 


where  a*   is  an  unbiased  estimate  of  a.   The  empirical  Bayes 
estimate  of  0,  commonly  called   the   mean   of   the  posterior 
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where  a*  is  an  unbiased  estimate  of  a.  The  empirical  Bayes 
estimate  of  9,  commonly  called  the  mean  of  the  posterior 
distribution  becomes 

5"  =  (..  fdl)   >»  +  '  (a'  t'd.)   >v.  i-i.a *   (3-5) 

This  is  the  mean  of  the  conditional  distribution  of  0  given 
x  [Ref .  8:p.  556] .  The  weights  of  each  of  the  prior  and 
sample  means  are  determined  by  the  prior  and  observed 
variances . 

As  a  measure  of  the  worth  of  this  estimator,  when  the 
number  of  x's  is  >  4  [Ref.  11] ,  the  risk,  or  sum  of  squared 
differences  between  the  true  unknown  parameter  and  the 
empirical  Bayes  estimate,  is  less  than  that  using  the 
observed  x's,  for  all  0i  .  Efron  and  Morris  [Ref.  9]  provide 
a  rigorous  proof. 

An  important  result  of  the  empirical  Bayes  estimator  is 
that  shrinkage  values  are  no  longer  constant,  but  depend 
upon  information  contained  in  each  cell.  Additionally,  the 
point  toward  which  the  cell  values  are  shrunk  is  the 
weighted  mean,  not  the  grand  mean  of  the  observations. 
Appendix  A  provides  the  derivation  of  the  empirical  Bayes 
estimator,  variable  shrinkage  factors  and  calculation  of 
the  prior  variance  used  in  this  study  as  they  pertain  to  the 
Marine  Corps  officer  attrition  rates. 


27 


IV .  VALIDATION  PROCEDURES  AND  RESULTS 

A.  GENERAL 

The  outputs  displayed  in  this  chapter  are  excerpts  from 
that  produced  by  the  Marine  Corps  Personnel  Inventory 
Attrition  Analysis  (MCIAAP)  (see  Appendix  B).  Initially, 
comparisons  will  be  made  with  the  output  displayed  by 
Robinson  to  identify  the  validity  of  the  Fortran  code  with 
the  APL  code  used  by  Robinson,  then  new  approaches  will  be 
displayed  in  an  attempt  to  improve  the  results.  The  Fortran 
code  in  the  MCIAAP  incorporates  most  of  the  procedures 
studied  by  Tucker  and  Robinson.  Specifically,  the  MLE , 
Transformed  Scale  Cell  Average  (TSCA)  and  the  James-Stein 
(JS)  estimators,  shrinkage  procedures  and  a  basic  inversion 
formula.  The  FTE  inversion  formula,  the  new  measures  of 
effectiveness  from  NPRDC ,  are  also  included.  For  the  MAD 
measure  of  effectiveness,  only  maximum  likelihood  and  the 
modified  empirical  Bayes  variables  will  be  presented,  since 
these  showed  the  most  promise  over  all.  All  other  variables 
are  displayed  by  the  MCIAAP  if  interest  is  warranted  by  the 
reader.  Chapter  II  describes  all  of  these  methods  in  some 
detail . 

B.  HISTORICAL  COMPARISONS 

Tables   4-9   were   developed   using  the  same  aggregation 
scheme  and  format  as   Robinson  [Ref .   2 : pp .  39-41],  although 
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net  all  of  the  variables  are  included  (MLE,  TSCA.  and  James 
Stein  are  the  only  variables  duplicated  in  this  paper) . 
Additionally,  Robinson  presented  the  number  of  cells,  as  k 
values,  used  for  evaluation  of  each  estimate  in  the 
aggregate  Figure  of  Merit  (FOM)  display.  These  k  values  are 
not  repeated  in  the  current  FOM  tables,  but  are  included,  as 
parenthetical  values,  and  discussed  in  the  new  MAD  measure 
of  effectiveness  tables.  Since  this  is  the  case  here,  the 
reader  is  reminded  that  the  apparent  competitive  performance 
of  the  MLE  in  the  following  cases  is  due,  in  part,  to  the 
possible  advantage  of  having  a  smaller  number  of  cells  used 
in  this  MOE ' s  computation. 

Since  the  computations  were  identical  for  both  studies 
for  these  tables,  the  output  was  expected  to  be  nearly 
identical  or  some  pattern  discernable  between  the  two.  In 
very  few  cases  was  this  found  to  be  true,  and  after 
verifying  that  the  algorithms  coded  both  in  the  APL  and 
Fortran  versions  were  correct,  we  have  no  reasonable 
explanation  for  the  discrepancies  of  the  programs.  It  is 
suspected  that  the  transfer  of  raw  data  into  an  APL 
workspace,  accomplished  by  Tucker  and  Robinson,  may  be  at 
variance  with  that  into  Fortran  files,  done  by  contractor. 

1.   Aviation  Aggregate  Case 

In  the  case  of  aviation  FOM ' s ,  displayed  in  Table  4, 
1981  values  for  LtCol ' s  in  the  transformed  scale  are 
identical  to  Robinson's  [Ref.  2:p.  39]  while  the  last  two 
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Table  4.     AVIATION  FIGURES  OF  MERIT 

1981 

1982 

1983 

t 

TRANSFORMED  FOM 

1st  Lt 

MLE 

4.0S63 

10.6590 

11.1338 

TSCA 

3.7420 

10.3636 

10.8711 

JS 

3.9694 

10.5671 

11.1626 

MODJS 

4.7352 

11.2615 

12.1206 

• 

EB 

4.0753 

10.7015 

11.3303 

MOD  EB 

3.1833 

5.0568 

6.2870 

LtCol 

MLE 

4.3664 

8.5033 

8.1506 

TSCA 

5.7768 

10.2994 

10.3353 

JS 

5.7369 

10.2468 

10.2972 

MOD  JS 

5.7208 

10.1865 

10.2756 

EB 

5.84S2 

10.2298 

10.3534 

MOD  EB 

4.3301 

3.9824 

2.7186 

ORIGINAL  FOM 

1st  Lt 

MLE 

21.0501 

43.1546 

53.0538 

TSCA 

56.4603 

93.7810 

98.4715 

TSCAM 

53.7S84 

90.5164 

95.3329 

JS 

58.0150 

95.2150 

100.1664 

JSM 

55.3844 

91.9406 

97.0163 

MODJS 

63.0345 

99.S419 

105.4897 

MOD  JSM 

60.6413 

96.5615 

102.3172 

EB 

5S.SS24 

96.2449 

101.24S2 

EBM 

56.2111 

92.9699 

98.0958 

MOD  EB 

19.7683 

38.4637 

49.7393 

MOD  EBM 

20.7987 

36.3673 

48.1659 

Lt  Col 

MLE 

33.9315 

59.4200 

22.1156 

TSCA 

38.4478 

57.0987 

69.0641 

TSCAM 

35.0390 

48.9374 

59.6957 

JS 

37.8248 

55.1913 

67.1696 

JSM 

33.2016 

46.9901 

56.7135 

MOD  JS 

38.0291 

51.5325 

63.3908 

MOD  JSM 

34.5466 

45.1554 

52.8285 

EB 

40.3421 

54.9389 

67.7091 

EBM 

36.7498 

46.9657 

57.8567 

MOD  EB 

40.2056 

70.5381 

60.5749 

MOD  EBM 

185.7476 

65.1553 

168.3976 
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years  are  smaller.  The  IstLt  FOM ' s  are  larger  for  all  three 
years.  The  modified  empirical  Bayes  (Mod  EB)  values  were 
best  overall  yet  still  may  be  too  large  (ideal  values  tend 
to  be  close  to  unity) .  It  is  interesting  to  note  that  in 
most  cases,  the  risk  values  tend  to  increase  with  future 
forecasts,  yet  the  Mod  EB  FOM  decreases  as  future  forecasts 
increase  with  time. 

In  the  original  scale,  MLE  tended  to  be  smaller  than 
Robinson's  for  all  years,  while  TSCA  and  James-Stein  (JS) 
were  all  larger.  The  FTE  inversion  worked  best  for  all 
variables  except  for  Mod  EB  in  the  first  and  third 
validation  years.  MLE  and  Mod  EB  worked  best  for  IstLt '3 
and  are  approximately  within  the  range  of  one  standard 
deviation  of  the  number  of  cells  in  the  computation.  For 
LtCol ' s ,  FTE  inverted  JS  was  best  with  MLE  a  competitive 
second,  yet  all  were  not  in  the  desired  range. 

In  the  MAD  measure,  Table  5,  results  are  mixed 
between  future  forecasts  (no  set  pattern)  and  variables  for 
the  two  ranks  (FTE  Mod  EB  was  best  for  IstLt ' s  and  MLE  was 
best  LtCol's).  The  mix  in  overage  and  underage  is  most 
desirable  in  the  first  year  for  all  cases,  yet  the  number  of 
cells  (k)  forecasted  as  overage  exceed  the  underage  forecast 
cell  numbers  in  most  cases.  This  is  cause  for  concern  and 
the  performance  of  these  estimates  in  this  aggregate  is 
suspect.  The  second  year  values  for  LtCol  estimates 
displays  a  preferred  ratio  of  k  values,  but  the  total  MOE 
value  may  still  be  too  large. 
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Table  5.     MEAN  ABSOLUTE  DEVIATION  (AVIATION) 

1981 

1982 

1983 

1st  Lt 

MLE 

UNDERAGE 

OVERAGE 

TOTAL 

0.1256(  3) 
2.3985(  8) 
2.5240(11) 

0.0439(  1) 
3.9267(10) 
3.9706(11) 

0.3069(  1) 
4.2522(10) 
4.5591(11) 

MOD  EB 

* 

UNDERAGE 

OVERAGE 

TOTAL 

0.1822(4) 
2.3625(  8) 
2.5447(12) 

0.04 10(  1) 
3.9040(11) 
3.9451(12) 

0.328I(  1) 
4.2351(11) 
4.5633(12) 

MOD  EBM 

UNDERAGE 

OVERAGE 
TOTAL 

0.2518(5) 
2.1654(  7) 
2.4172(12) 

0.0580(  2) 
3.6427(10) 
3.7007(12) 

0.3420(  1) 
3.9677(11) 
4.3098(12) 

LtCol 

MLE 

UNDERAGE 

OVERAGE 

TOTAL 

0.6005(  5) 

0.7925(  8) 
1.3930(13) 

1.5139(10) 
0.3928(  3) 
1.9067(13) 

0.0396(  1) 

1.3467(12) 
1.3863(13) 

MOD  EB 

UNDERAGE 
OVERAGE 

TOTAL 

0.6053(  6) 

1.3083(10) 
1.9136(16) 

1.4717(  9) 
1.0970(  8) 
2.56S7(17) 

0.0841  (  2) 
2.2331(14) 
2.3173(16) 

MOD  EBM 

UNDERAGE 

OVERAGE 

TOTAL 

0.6890(  6) 
1.0120(10) 
1.7010(16) 

1.5366(10) 
0.71 79(  7) 
2.2546(17) 

0.1 167(  2) 
1.7501(14) 
1.8669(16) 

2.   Combat  Support  Aggregate  Case 

Table  6  displays  the  combat  support  aggregate  output 
from  the  program.  In  the  transformed  scale  comparison  with 
Robinson  [Ref.  2:p.  401,  the  values  for  the  two  ranks  take 
opposite  directions.   Where  the  IstLt  FOM  values  are  all 
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Table  6.     COMBAT  SUPPORT  FIGURES  OF  MERIT 

1981 

1982 

1983 

TRANSFORMED  FOM 

l 

1st  Lt 

MLE 

2.1877 

2.9028 

2.2447 

TSCA 

1.6329 

2.4580 

1.8696 

JS 

1.6925 

2.4816 

1.8987 

MOD  JS 

2.3684 

3.0520 

2.5356 

EB 

2.3323 

3.1991 

2.6S68 

MOD  EB 

2.2293 

2.0048 

'      2.1531 

LtCol 

MLE 

1.3103 

1.6978 

2.3108 

TSCA 

0.7502 

1.4070 

1.7540 

JS 

0.6421 

1.2345 

1.5429 

MOD  JS 

0.7009 

1.0929 

1.3082 

EB 

0.7736 

1.25S9 

1.5799 

MOD  EB 

1.S691 

2.1236 
ORIGINAL  FOM 

4.9017 

1st  Lt 

MLE 

111.4260 

112.1917 

70.2767 

TSCA 

73.5008 

90.0523 

73.93S0 

TSCAM 

79.0947 

93.8244 

80.9764 

JS 

76.9945 

89.8064 

71.5890 

JSM 

429.8503 

84.4175 

65.2639 

MOD  JS 

143.1700 

113.0083 

90.6783 

MOD  JSM 

168.9385 

114.7083 

92.6599 

EB 

102.9525 

119.0900 

96.3373 

EBM 

108.3578 

111.2701 

86.8138 

MOD  EB 

6S.7251 

63.9054 

48.2719 

MOD  EBM 

176.3961 

59.5760 

43.0037 

LtCol 

MLE 

36.0152 

66.0261 

42.4097 

TSCA 

61.7769 

52.0346 

99.7917 

TSCAM 

65.9535 

49.6307 

107.1016 

JS 

61.0511 

49.2280 

96.7861 

JSM 

66.3509 

47.8490 

118.2564 

MOD  JS 

66.5470 

49.9366 

99.2566 

MOD  JSM 

81.4722 

50.2637 

103.8460 

EB 

69.7508 

55.9812 

104.5311 

EBM 

71.5728 

50.0281 

109.2012 

MOD  EB 

62.6316 

50.3125 

96.8693 

MOD  EBM 

83.1337 

4S.6576 

127.1541 
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larger  (except  for  MLE  in  1981) , the  LtCol  values  ar= 
smaller.  Overall,  all  estimates  are  competitive  with  each 
other  over  all  validation  years  with  TSCA  performing  best 
for  the  IstLt's  and  James-Stein  outperforming  the  rest  in 
the  LtCol's.  Note  here  also  that  the  Mod  EB  estimate  did 
not  fair  well  in  either  case,  and  was  the  worst  estimate  in 
risk  for  the  LtCol's. 

In  the  original  scale  ,  a  different  and  more  random 
effect  is  noticed  between  Table  6  and  the  results  of 
P.obinson.  Risk  values  tended  to  fluctuate  with  forecast 
lead  time  between  estimates  for  IstLt's  where  the  LtCcl 
values  were  more  stable,  yet  were  quite  different  from  the 
transformed  scale  pattern.  In  this  case,  Mod  EB  performed 
well  for  IstLt's,  and  was  the  only  one  within  the  desired 
range  of  standard  deviations.  MLE  for  the  LtCol's  was  also 
the  only  estimate  desirable  for  this  measure.  The  Miller 
inverse  showed  poorly  in  this  aggregate  for  all  estimation 
methods,  in  every  case  exceeding  it's  counterpart. 

In  the  MAD  measure  of  effectiveness,  Table  7,  MLE 
performed  the  best  in  both  ranks.  The  Mod  EB  effect  is 
competitive  in  the  IstLt's  case,  and  all  variables  in  the 
LtCol  case  are  closely  grouped.  The  number  of  cells  are 
heavily  weighted  towards  an  overage  forecast  in  all  three 
validation  years,  throughout  this  aggregate.  Once  again, 
this  is  not  a  desirable  quality  for  a  MOE . 
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Table  7.     MEAN  ABSOLUTE  DEVIATION  (COMBAT  SUPPORT) 


19S1 

1982 

1983 

1st  Lt 

MLE 

UNDERAGE 

OVERAGE 
TOTAL 

0.3737(14) 
0.6731(19) 
1.0468(33) 

0.1617(8) 

0.9574(25) 
1.1191(33) 

0.1781(8)      • 

0.8835(25) 

1.0616(33) 

MOD  EB 

UNDERAGE 

OVERAGE 
TOTAL 

0.4028(15) 
0.7311(40) 
1.1339(55) 

0.1275(9) 
0.9972(45) 
1.1246(54) 

0.1427(9) 
0.9422(38) 
1.0849(47) 

MOD  EBM 

UNDERAGE 

OVERAGE 
TOTAL 

0.4592(16) 
0.5930(39) 
1.0522(55) 

0.1603(9) 
0.S3 18(45) 
0.9921(54) 

0.1856(10) 

0.7575(37) 
0.9431(47) 

LtCol 

MLE 

UNDERAGE 
OVERAGE 

TOTAL 

0.1360(  9) 
0.1510(14) 
0.2S70(23) 

0.1688(  7) 
0.2591(16) 
0.4280(23) 

0.054 1(  5) 
0.3415(18) 
0.3957(23) 

MOD  EB 

UNDERAGE 

OVERAGE 

TOTAL 

0.1638(11) 
0.2232(29) 
0.3S70(40) 

0.1721(  7) 
0.3606(35) 
0.5327(42) 

0.1 162(  6) 

0.4548(29) 
0.5710(35) 

MOD  EBM 

UNDERAGE 

OVERAGE 

TOTAL 

0.2044(11) 
0.1263(29) 
0.3307(40) 

0.1901(  7) 
0.2240(35) 
0.4142(42) 

0.140S(  6) 
0.2589(29) 
0.3996(35) 
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3.   Ground  Combat  Aggregate  Case 

In  the  final  aggregate  studied  by  Rcbinson  [Ref. 
2:p.  41] ,  that  of  ground  combat  MOS ' s ,  the  results  again,  as 
displayed  in  Table  8,  were  mixed  between  computational 
methods.  In  the  original  scale,  value  differences  tended  to 
be  larger  than  Robinson's  in  all  cases  for  IstLt's,  showing 
a  proportional  growth  over  time.  Conversely,  LtCcl  values 
were  smaller  in  the  first  two  validation  years  and  slightly 
larger  in  the  final  year.  TSCA  looks  good  for  the  IstLt's 
with  the  others  being  competitive,  yet  the  size  of  the 
values  may  net  be  desirable.  Mod  EB  turned  out  the  best  for 
LtCol ' s ,  but  all  other  estimates  are  competitive  and  in  the 
desired  range. 

On  the  original  scale,  all  values  for  the  estimates 
for  IstLt's  are  too  large,  with  Mod  EB  performing  well  in 
that  regard.  LtCol  values  are  all  larger  than  Robinson's, 
but  the  estimates  are  all  within  range  of  one  standard 
deviation  and,  therefore,  competitive.  MLE ,  once  again, 
shows  well  for  this  rank.  In  opposing  fashion  to  that  of 
the  combat  support  aggregate,  the  FTE  inversion  outperformed 
all  Basic  inversions  for  all  variables. 

The  MAD  measure  of  effectiveness,  Table  9  confirms 
the  results  found  above  in  that  Mod  EB  was  best  for  IstLt's 
and  MLE  did  well  for  LtCol ' s .  The  values  for  IstLt  were 
found  to  be  too  large  overall,  while  all  estimates  for  LtCol 
were  competitive   and  in  the  desired  range.   The  mix  between 
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Table  8.     GROUND  COMBAT  FIGURES  OF  MERIT 

1981 

19S2 

1001 

iVoj 

« 

TRANSFORMED  FOM 

1st  Lt 

MLE 

2.5764 

4.4136 

6.8744 

TSCA 

1.9741 

3.8072 

6.H52 

JS 

2.1737 

3.9738 

6.3970 

MOD  JS 

3.0859 

4.7919 

7.5619 

EB 

2.5558 

4.3626 

6.9607 

MOD  EB 

3.1538 

3.5202 

4.7069 

LtCol 

MLE 

1.2126 

1.9674 

3.1084 

TSCA 

1.5715 

2.3553 

3.5566 

JS 

1.5571 

2.2719 

3.4365 

MOD  JS 

1.6227 

2.1317 

3.1S83 

EB 

1.7087 

2.2464 

3. j  j-o 

MOD  EB 

1.1964 

I 

2.5830 
ORIGINAL  FOM 

2.5029 

1st  Lt 

MLE 

87.3362 

119.2637 

143.9128 

TSCA 

87.9712 

116.6767 

145.6236 

TSCAM 

83.7648 

111.3995 

142.0902 

JS 

95.4717 

122.9425 

152.7354 

JSM 

92.0238 

117.5264 

149.3965 

MOD  JS 

124.6490 

147.6611 

180.5502 

MOD  JSM 

134.9472 

142.5461 

178.6173 

EB 

109.6349 

136.6321 

167.3234 

EBM 

105.5708 

131.0400 

163.S189 

MOD  EB 

64.3647 

81.2071 

101.9511 

MOD  EBM 

63.6623 

77.S159 

100.7150 

LtCol 

MLE 

35.0538 

45.8172 

53.0111 

TSCA 

55.9649 

198.8659 

123.1306 

TSCAM 

43.5234 

201.9953 

122.8648 

JS 

56.3346 

198.9037 

122.2786 

JSM 

43.9957 

201.3299 

115.7143 

MOD  JS 

60.0730 

202.2959 

122.9959 

MOD  JSM 

52.3940 

210.7551 

121.4254 

EB 

62.7441 

205.4249 

126.6544 

EBM 

50.4049 

209.3853 

125.7642 

MOD  EB 

48.0871 

197.7316 

111.9375 

MOD  EBM 

3S.1231 

323.5654 

242.4398 
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over-estimates  and  under-estimates ,  i.e.,  the  number  of 
cells  forecasted  over  or  under  the  actuals,  is  not  a 
desirable  one,  yet  the  same  mix  holds  true  for  all 
validation  years.  This  also  was  the  finding  of  both 
preceding  aggregates. 


Table  9.     MEAN  ABSOLUTE  DEVIATION  (GROUND  COMBAT) 


1981 

1982 

1983 

1st  Lt 

MLE 

UNDERAGE 

OVERAGE 

TOTAL 

0.1525(5) 
1.9322(30) 

2.0847(35) 

0.1855(10) 

2.5979(25) 
2.7834(35) 

0.2466(10) 
3.1691(25) 

3.4157(35) 

MOD  EB 

UNDERAGE 

OVERAGE 

TOTAL 

0.2340(10) 

1.7743(33) 
2.0083(43) 

0.1323(  8) 
2.3086(32) 
2.4409(40) 

0.2220(  8) 
2.8449(30) 
3.0669(38) 

MOD  EBM 

UNDERAGE 

OVERAGE 
TOTAL 

0.2618(10) 

1.6492(33) 
1.9111(43) 

0.1705(10) 
2.173S(30) 
2.3443(40) 

0.2617(10) 

2.6935(28) 
2.9552(38) 

LtCol 

MLE 

UNDERAGE 

OVERAGE 

TOTAL 

0.2313(10) 
0.4182(25) 
0.6495(35) 

0.3034(11) 
0.4202(24) 
0.7236(35) 

0.136S(10) 
0.6453(25) 
0.7821(35) 

MOD  EB 

UNDERAGE 

OVERAGE 

TOTAL 

0.2315(11) 
0.6392(34) 
0.S706(45) 

0.3299(14) 
0.6922(32) 
1.0221(46) 

0.1426(12) 
0.9523(31) 
1.0950(43) 

MOD  EBM 

UNDERAGE 

OVERAGE 

TOTAL 

0.2650(12) 
0.4071(33) 
0.6720(45) 

0.3790(15) 

0.4427(31) 
O.S217(46) 

0.1864(12) 
0.6535(31) 
0.8399(43) 
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4.   Test  Aggregate  Case 

In  an  effort  to  experiment  with  aggregate  groupings, 
another  case  was  built  upon  which  all  estimation  methods 
were  applied,  and  the  results  compared  with  previously 
studied  groups.  The  grade  was  held  fixed,  in  this  case  the 
grade  of  Captain  was  used,  with  a  LCS  in  the  eighth  year. 
Several  MOS  '  s  were  chosen  so  that  to  insure  an  aggregate 
with  non-zero  inventories.  In  the  transformed  scale,  Table 
10,  all  estimates  have  small  and  competitive  values,  with 
Mod  EB  and  MLE  doing  well.  In  the  original  scale,  once 
again,  all  estimates  are  well  within  range,  but  Mod  EB  and 
MLE  were  the  standouts.  The  FTE  inversion  option  did  not 
fair  well  with  this  aggregate.  All  estimates  in  the  MAD 
measure  of  effectiveness,  Table  11,  are  very  competitive, 
and  the  first  validation  year  mix  (overage  and  underage)  may 
be  the  best  yet.  It  is  interesting  to  note  that  the  cell 
split  in  the  first  validation  year  between  overage  and 
underage  forecasts  is  contrary  to  that  of  the  historical 
displays  above.  Although  the  total  number  of  cells  used 
(k=6)  is  much  smaller  than  the  previous  table,  it  is  felt 
that  the  aggregate  scheme  was  the  major  factor  in  this 
finding,  and  the  results  are  promising. 


39 


Table  10.     TEST  AGGREGATE  FIGURES  OF  MERIT 

1981 

1982 

1983 

TRANSFORMED  FOM 

Captain 

• 

■ 

MLE 

0.9603 

2.3463 

2.0842 

TSCA 

1.4151 

3.3033 

2.9152 

JS 

1.3753 

3.3969 

2.8725 

MOD  JS 

1.3943 

3.S165 

2.8829 

EB 

1.6S40 

4.2698 

3.2957 

MOD  EB 

0.8249 

1.7132 
ORIGINAL  FOM 

1.29S0 

Captain 

MLE 

3.0402 

4.8104 

4.9399 

TSCA 

4.6161 

10.2779 

7.5714 

TSCAM 

8.1467 

7.6241 

8.8347 

JS 

4.3970 

10.3000 

7.4726 

JSM 

7.3199 

7.6666 

S.3S60 

MOD  JS 

4.2990 

10.7640 

7.6371 

MOD  JSM 

6.0079 

8.3070 

7.9259 

EB 

5.3660 

12.0337 

S.8679 

EBM 

6.7181 

9.3915 

8.7959 

MOD  EB 

1.4463 

6.0323 

3.9066 

MOD  EBM 

4.6465 

4.3354 

5.9175 
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Table  11.     MEAN  ABSOLUTE  DEVIATION  (TEST  CASE) 

1981 

1982 

19S3 

Captain 

, 

• 

MLE 

UNDERAGE 

OVERAGE 

TOTAL 

0.2625(  3) 
0.1351(  3) 
0.3976(  6) 

0.0463(  1) 
0.7374(  5) 
0.7837(  6) 

0.2097(  2) 
0.4391(  4) 
0.6488(  6) 

MOD  EB 

UNDERAGE 

OVERAGE 

TOTAL 

0.1342(4) 
0.2575(  2) 
0.39 18(  6) 

0.0000(  0) 
0.9601(  6) 
0.9601(  6) 

0.0965(  1) 
0.6306(  5) 
0.7270(  6) 

MOD  EBM 

- 

UNDERAGE 

OVERAGE 

TOTAL 

0.3712(  4) 
0.0852(  2) 
0.4563(  6) 

0.0809(  1) 
0.6148(  5) 
0.6957(  6) 

0.2352(  2) 
0.3025(  4) 
0.5377(  6) 

.   * 
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V.   CONCLUSIONS  AND  RECOMMENDATIONS 

A.   RESULTS 

Omitting  the  fact  of  differences  in  computational 
methods  between  this  paper  and  the  Robinson  thesis,  it  is 
clear  that  the  results  are  still  inconclusive  at  this  time. 
In  the  historical  comparisons,  the  estimates  studied  varied 
in  their  performance  from  one  aggregate  to  another.  The 
different  risk  calculations,  in  the  transformed  and  original 
scales,  as  well  as  the  different  measures  of  effectiveness 
showed  little  agreement  in  estimation  methods.  Long  range 
estimation  still  remains  a  problem  since  with  time,  the 
estimate  risks  increase  significantly.  Much  emphasis  was 
put  on  the  modified  empirical  Bayes  estimator  in  this  paper 
and  those  specific  results  look  promising  under  possibly 
different  and  more  stable  aggregate  conditions. 

One  area  of  mention  here  concerning  the  results  has  to 
do  with  the  composition  of  inventories  used  in  the  selected 
aggregates.  It  was  noticed  that  the  range  of  inventory 
(smallest  to  largest  value  in  any  particular  MOS )  was  quite 
large,  and  may  be  affecting  the  results.  In  fact,  Carter 
and  Rolph  [Ref.  4:p.  883]  managed  this  point  by  further 
disaggregating  into  levels  of  activity  as  it  concerned  their 
work  with  fire  alarm  estimation. 
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B.  CONCLUSIONS 

The  results  of  the  previous  feasibility  studies  and  this 
paper  leads  to  several  conclusions.  The  chosen  aggregate 
still  seems  to  be  the  key  element  in  the  work  of  attrition 
estimation.  It  does  seem  clear  that  the  aggregation  scheme 
should  be  focused  on  fixed  ranks,  in  narrow  LOS  ranges  and 
multiple  MOS's.  This  conclusion  is  based  upon  the  relative 
success  of  the  test  aggregate  utilized  in  this  work,  and 
upon  similar  success  by  Carter  and  Rolph.  Additionally,  the 
partitioning  of  inventory  levels  in  the  aggregate  should 
stabilize  the  estimation  method  performance.  All  the 
estimation  method  studied  to  date  have  their  strengths,  but 
still,  no  one  method  of  estimation  excels  over  the  others. 

C.  RECOMMENDATIONS 

No  estimation  method  studied  in  this  paper  is 
recommended  for  implementation  at  this  time.  In  addition, 
no  method  is  recommended  for  removal  from  consideration 
since  basic  aggregation  problems  are  felt  to  be  the  largest 
contributory  factor  to  the  mixed  results  found  here.  These 
recommendations  are  presented  for  further  study. 

1.  The   aggregation   method   offered   by   Captain  Larsen 

[Ref.   3]    should   be    examined   further   with   its 
relevance  to  the  scheme  in  the  test  case  used  here. 

2.  Utilization  of  the  detailed  data  provided  by  NPRDC  as 
soon  as  possible  in  order  to  expand  the  available 
information  and  open  new  avenues  of  aggregation. 

3.  Re-examine  the  inventory  mean  value  used  in  the  FTE 
inverse  to  determine  its  future  role  in  the 
estimation  process. 
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APPENDIX  A.   ESTIMATION  ALGORITHMS 

A.   NOTATION 

The  following  notation  is  a  generalization  of  that  used 
by  Robinson  [Ref.  2:p.  71]  as  it  pertains  to  the  aggregation 
scheme  used  in  this  paper.  The  indexing  system  has  been 
changed  to  reflect  the  use  of  officer  grades  to  identify  the 
cells  . 

i  =  index  of  MOS  cells  in  the  aggregate. 

j  =  index  of  LOS  cells  in  the  aggregate. 

k  =  index  of  Grade  cells  in  the   aggregate. 

t  =  index  of  time  periods  in  the  aggregate. 

inv ( i , j , k r t )  =  inventory  with  MOS  i,  LOS  j  and  Grade  k 

during  year  t. 
y(i,j,k,t)  =  number  of  attritions  in  cell  (i,j,k)  during 

year  t. 
n ( i , j , k , t )  =  maximum  i y ( i , j , k , t ) , . 5 [inv ( i , j , k , t )  + 

inv(i, j ,  k, t  +  1)  ]  I 
D  =  matrix  identifying  cells  with  no  inventory  over  all 

estimatior  years. 
d(i,j,k)  =  0,  if  cell  is  a  structural  zero. 
d(i,j,k)  =  1,  if  cell  is  not  a  structural  zero. 
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B.   FREEMAN-TUKEY  EXACT  INVERSION  ALGORITHM 

This  algorithm  provides  a  formula  for  the  inverse  of  the 
Freeman-Tukey  double  arcsine  transformation  to  express  the 
means  of  the  transformed  values  into  proportions  on  the 
original  scale. 

STEP  1.    Transform  the  data  using  the  Freeman-Tukey 
double  arcsine  transformation. 

x(i,j,k,t)  =  0.5[n(irj,krt)+0.5p'2  (sin*  ^  (^  ( I (  j  '  k  !  t ) +1 )  "1] 

+  sin-  E?*??1'^*;^  "^  I  <A.l) 

(n  (1 ,  3 , k, t) +1) 

STEP  2.   Calculate   the   average   cell  inventory  and 
transformed  values  over  all  T. 


x(i,j,k)  =  (1/T)S  x(i,j,k),  (A. 2) 

t 


inv(i,j,k)  =  (1/T) S  inv (i , j ,k) .       (A. 3) 

t 


STEP  3.  Check  for  inversion  values  less  than  zero 
or  greater  than  one.  The  algorithm  is  terminated  at  this 
point  if  these  conditions  are  met. 


tm  =  x(i, j,k) /(inv(i, j ,k)+0.5M /2  +  n/2,       (A. 4) 


tmi  =  sin-1 (1/ (inv(i, j,k)  +  l))1/2,  (A. 5) 

If  tm  <  tm! ,         FTE  =  0.0 , 
IF  tm  >  (n-tmi ) ,     FTE  =  1.0, 
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STEP  4.   Compute  the  sgn  function. 

If  cos(tm)  >  O.C     sgn  =  1.0, 
If  cos(tm)  <  0.0     sgn  =  -1.0 


STEP  5.   Computer  the  inverted  value  of  the 
transformed  variable 

FTE  =  0.5(l-sgnll-[sin(tm)  +  sin (tm)  -  1/sin ( tm) j 2  ji/a  , 

mv(i,:,k) 


C.   TRANSFORM  SCALE  VARIANCE  ALGORITHM 

This  algorithm  produces  the  variance  of  the  transformed 
data  within  the  unstable  variance  range  where  the  assumption 
of  normality  does  not  hold.  As  discussed  in  Chapter  II, 
this  range  of  instability  was  determined  to  be  between  1.001 
and  2.2.  The  method  for  computing  the  variance  within  this 
range  is  presented  here. 

STEP  1.  Compute  the  transformed  scale  values  to  be 
used  for  the  variance.  Since  x(i,j,k)  was  computed  in  the 
FTE  inversion  algorithm  of  this  appendix,  the  formulae  will 
not  be  repeated. 


x(t)  =  z:  [0.5  +  inv(i,  j  ,  k,  t)  I1  >'  =  ,  (A.  7) 

ijk 


z(i,j,k)  =  x(i,j,k)  +^r^-(n/2).  (A. 3) 
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STEP  2.  Compute  the  variance  of  transformed  cell 
values.  The  formula  is  derived  in  Chapter  II  as  a  linear 
regression  equation.  The  range  of  the  variance  is  limited 
to  0.05  <  var(i,j,k)  <,    1.0. 

var(i,j,k)  =  az  (i  ,  j  ,  k) h  l   (z  (i  ,  j  ,  k)  -  IP2, 

1.001  <z(i,j,k)<;  2.2  (A. 9) 

where  a  =  1.6335,  bl  =  -0.3934,  and  b2  =  0.8991. 

D.   EMPIRICAL  BAYES  SHRINKAGE  ALGORITHM 

This   algorithm   computes   the  empirical  Bayes  shrinkage 
values  and  estimator  of  attrition  rates.   Prior   to  entering 
this  algorithm,  the  value  for  a  has  been  set  to  0.0. 
STEP  1.   Initialization. 

ai  =  a.  (A. 10) 

STEP  2.  Define  the  a  and  T  parameters  in  the 
empirical  Bayes  context. 


a  =  H    l/(a  +  var  (i  ,  j  ,  k)  )  , 
ijk 

where  var(i,j,k)  is  defined  by  formula  (A. 9) 


r=    H    a/    11    1/  (a    +    var  ( i  ,  j  ,  k }  )  .  ( A .  12 ) 

ijk  ijk 
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STEP  3.   Compute  the  empirical  Bayes  prior  variano 


k-1-   21  a[x(i, j ,  k) -z]2 

a  =  a  -  ijk ,  (A. 13) 

21  [a  (x(i,  j  ,  k)  -  z]  - 
ijk 


where    z    =    21    ^(i,j,k)     and    k    =    11    d(i,j,k) 
ijk  ijk 


STEP  4.   Check   for   the   following   conditions   and 
branch  accordingly. 

If,  a  <  0.0  set  a  =  0.0  and  go  to  step  5, 

If,  |a  -  ai I  >  .0001     return  to  Step  1. 


STEP  5.   Compute  the  empirical  Bayes  estimator 


::E  B  (  i  ,  j  ,  k )  =  [— 7-: — — 

a+var(i,j, 


k 


]x(i,j,k)  +        (A. 14) 


r    var(i,,j,k)  ,- 
La+var (i, j,k) J 


STEP  6.   Compute  the  empirical  Bayes  shrinkage  set. 

shrEB(i,j,k)  =  var (i, j , k) / (a  +  var(i,j,k)}.       (A. 15) 

E.   MODIFIED  EMPIRICAL  BAYES  SHRINKAGE  ALGORITHM 

This  algorithm   differs  only   slightly  from   that  of  th- 
empirical  Bayes  shrinkage  algorithm  listed  above,  and  for 
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that   reason,   only   the   pertinent   steps  will  be  displayed 
here  . 

Consider  a  particular  cell  and  its  transformed  value 
x(t)  for  t=l,...,T  in  the  estimation  set.  To  form  a 
weighted  time  average,  set 

xt(t)  =  x(t)/[0.5  +  inv(t)]1/2  (A. 16) 

This  leads  to  a  modification  of  step  1  in  the  following 
form. 

STEP  1.   Compute   the    weighted   transformed   scale 
average . 

xt  =  i  S  xt(t)  (A. 17) 

t 

Since  the  weights  will  modify  the   variance  of   the  weighted 
transformed  scale  values. 

var[xt(t)]  =  var (i, j ,k) / (0.5  +  inv(t))         (A. 13) 

var[xt(t)]  =  (h2     ZI  var[xt(t)]  =  vt/T         (A. 19) 

t 

The  extreme  right  of  formula  A. 19  serves  to  define  vt , 
the  approximation  for  the  weighted  variance  in  the 
validation  set  for  that  cell. 

The  values  of  xt  and  var(xt)  are  carried  throughout  this 
version  of  empirical  Bayes  estimation  to  form  the  values 
XTEB  for  the  empirical  Bayes  attrition  rate.  For  the 
purposes  of  inversion  to  the  original  scale,  the  usual  Basic 
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formula  is  modified  since  the  weighting   factor,  0.5   +  inv, 
is  already   incorporated  into   the  transform  value.   The  new 
Basic  inversion  for  this  version  is  simply  as  follows: 

p  =  1/2[1  +  sin(XTEB)]  (A. 20) 
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APPENDIX  B.     MARINE  CORPS  PERSONNEL  INVENTORY 
ATTRITION  ANALYSIS  PROGRAM 

PROGRAM  MARCOR 

* 

* 
* 

* 
Vr 

ft 

* 


MARINE  CORPS   PERSONNEL  INVENTORY  ANALYSIS   PROGRAM 

PURPOSE:    TO  PROVIDE  FUTURE   PERSONNEL  LOSS  ESTIMATION  OF  THE 
MARINE  CORPS  OFFICER  STRUCTURE  BY  SEVERAL  DIFFERENT 
ESTIMATION  SCHEMES. 
PROGRAMMERS:    LUIS  URIBE,    INDEPENDENT  CONTRACTOR 

CAPT  C  R  DICKINSON,    USMC 
DESCRIPTION:    WRITTEN  IN  FORTRAN  77  FOR  THE   IBM  3033  MAINFRAME 
COMPUTER  RESIDENT  AT  NPS,   THIS   PROGRAM   INCORPORATES 
METHODS  EXPLORED  BY  MAJOR  D.  D.    TUCKER    ,    MAJOR  J.  R. 
ROBINSON  AND  CAPTAIN  DICKINSON.    ESTIMATION  SCHEMES 
INCLUDED  ARE: 

MAXIMUM  LIKELIHOOD  ESTIMATION   (MLE) 
TRANSFORMED  SCALE  CELL  AVERAGE   (TSCA) 
JAMES-STEIN 

LIMITED  TRANSLATION  JAMES -STEIN 
EMPIRICAL  BAYES   ESTIMATION 
RISK  ANALYSIS  WAS   PERFORMED   IN  BOTH  THE  ORIGINAL  AND 
TRANSFORMED  SCALES.    ON  THE  TRANSFORMED  SCALE, 
ACHIEVED  BY  USING  THE  FREEMAN-TUKEY  DOUBLE  ARCSINE 
TRANSFORMATION,   THE  RISK  WAS  DEFINED  BY  THE  AVERAGE 
SQUARED  DEVIATION  OF  ACTUAL  AND  ESTIMATED  VALUES. 
FOR  THE   ORIGINAL  SCALE,    THE  TRANSFORMED  VALUES  WERE 
FIRST   INVERTED  USING  TWO  TECHNIQUES.    AN  AD-HOC 
FORMULA  USED  BY  TUCKER  AND  ROBINSON  AND  AN  EXACT 
INVERSE  FORMULA  FROM  JOHN  J  MILLER.    TWO  MEASURES  OF 
EFFECTIVENESS,    CHI -SQUARE   GOODNESS   OF  FIT  AND  MEAN 
ABSOLUTE  DEVIATION  (MAD),   WERE  THEN  EVALUATED. 
RESULTS  ARE   PRINTED,   WITH  ALL  ARRAYS  COMPUTED,    AT  THE 
COMPLETION  OF  THE  PROGRAM. 
VARIABLE  DESCRIPTION: 

THE  FOLLOWING  LIST  OF  VARIABLES,   WITH  A  BRIEF 
DESCRIPTION,    FORM  THE   BASIS  FOR  THE  VARIABLE  NAMES 
FOR  THE  ENTIRE  PROGRAM. 

MLE   -  MAXIMUM  LIKELIHOOD  ESTIMATOR 
TSCA   -  TRANSFORMED   SCALE  CELL  AVERAGE 
J   -   JAMES-STEIN  ESTIMATOR 
Jl   -   JAMES-STEIN  ESTIMATOR  MODIFIED  BY 
A  VARIATION  ON  SHRINKAGE 
JLT   -   LIMITED  TRANSLATION  JAMES-STEIN 

EB    -   EMPIRICAL  BAYES 
TEB   -  EMPIRICAL  BAYES  MODIFIED  BY  A 
VARIATION  ON 
TRANSFORMED  SPACE  VARIABLES 

THE  FOLLOWING  ARE  UNIQUE  VARIABLES   IN  THE  TRANSFORMED 
SPACE: 

X  -  VALUE  OF  THE  FREEMAN-TUKEY 


MAR00010 
MAR00020 
MAR00030 
MAR00040 
MAR00050 
MAR00060 
MAR00070 
MAR00080 
MAR00090 
MAR00100 
MAR00110 
MAR00120 
MAR00130 
MAR00140 
MAR00150 
MAR00160 
MAR00170 
MAR00180 
MAR00190 
MAR00200 
MAR00210 
MAR00220 
MAR00230 
MAR00240 
MAR00250 
MAR00260 
MAR00270 
MAR00280 
MAR00290 
MAR00300 
MAR00310 
MAR00320 
MAR00330 
MAR00340 
MAR00350 
MAR00360 
MAR00370 
MAR00380 
MAR00390 
MAR00400 
MAR00410 
MAR00420 
MAR00430 
MAR00440 
MAR00450 
MAR00460 
MAR00470 
MAR00480 
MAR00490 
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TRANSFORMATION 
XB  -  (X-BAR)  AVERAGE  TRANSFORM  VALUE 
OVER  ESTIMATION  YEARS 
XBB  -  GRAND  MEAN  OF  TRANSFORMED  VALUES 
SSE  -  SUM  OF  SQUARED  DIFFERENCES,  ERROR 

(X  AND  X-BAR) 
SST  -  SUM  OF  SQUARED  DIFFERENCES,  TOTAL 

(X  AND  XBB) 
SSB  -  SUM  OF  SQUARED  DIFFERENCES,  BETWEEN 
(SST-SSE) 
SHR  -  JAMES-STEIN  SHRINKAGE  VALUES 
SHR1  -  MODIFIED  JAMES-STEIN  SHRINKAGE 
VALUES 
SHREB  -  EMPIRICAL  BAYES  SHRINKAGE  VALUES 
SHRTEB  -  MODIFIED  EMPIRICAL  BAYES  SHRINKAGE 
VALUES 
THE  FOLLOWING  ARE  MODIFIED  BASIC  VARIABLES  IN  THE 
TRANSFORMED  SPACE.  THIS  IS  DONE  BY  PLACING  AN  'X'  FOR 
A  PURE  TRANSFORMED  VARIABLE,  OR  AN  'R'  FOR  A  RISK 
VARIABLE  IN  FRONT  OF  THE  BASIC  VARIABLE. 
MLE  RML  (MLE) 

TSCA  RTSCA 

XJ  RJ 

XJ1  RJ1 

XJLT  RJLT 

XEB  RXEB 

XTEB  RTEB 

ORIGINAL  SPACE  VARIABLES 

THE  FOLLOWING  ARE  MODIFIED  BASIC  VARIABLES  IN  THE 
ORIGINAL  SPACE.  SIMILAR  TO  WHAT  IS  DONE  IN  THE 
TRANSFORMED  SPACE.  A  'P'  SIGNIFIES  A  PURE  INVERTED 
VARIABLE,  AND  AN  fRR'  REPRESENTS  A  RISK  VARIABLE.  AN 
'M'  AT  THE  END  OF  A  VARIABLE  NAME  IDENTIFIES  A  MILLER 
INVERSE  VARIABLE. 

PMLE  RRML  (MLE) 

PJ  RRJ 

PJM  RRPJM 

PJ1  RRJ1 

PJ1M  RRJ1M 

PJLT  RRJLT 

PEB  RRXEB 

PEBM  RRXEBM 

PTEB  RRRTEB 

PTEBM         RRTEBM 
THE  FOLLOWING  LIST  IS  OF  PREFIXES  USED  TO  DESCRIBE 
THE  MEAN  ABSOLUTE  DEVIATION  VARIABLES: 
MO  -  MEAN  DEVIATION  'OVERAGE* 
MU  -  MEAN  DEVIATION  'UNDERAGE' 
MD  -  MEAN  'TOTAL  DEVIATION 
KVO  -  NUMBER  OF  'OVERAGE'  ESTIMATES  IN  THE 

AGGREGATE 
KVU  -  NUMBER  OF  'UNDERAGE'  ESTIMATES  IN  THE 

AGGREGATE 
KVD  -  NUMBER  OF  'TOTAL'  ESTIMATION  DEVIATIONS 
IN  THE  AGGREGATE 
THESE  PREFIXES  ARE  APPENDED  TO  THE  BASIC  VARIABLES  TO 


* 
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■ft 
k 
■ft 
■ft 
k 
■ft 
k 
■ft 
it 
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k 
■k 
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k 
k 
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FORM  THE  MAD  VARIABLES  AS  FOLLOWS 

MOTSCA  MUTSCA  MDTSCA  KVOTSA  KVUTSA  KVDTSA  (TSCA) 

MOTSAM  MUTSAM  MDTSAM  KVOTSM  KVUTSM  KVDTSM  (TSCAM) 

MOPMLE  MUPMLE  MDPMLE  KVOPML  KVUPML  KVDPML  (PMLE) 

MOP J   MUPJ   MDPJ   KVOPJ  KVUPJ  KVDPJ  (PJ) 

MOPJM  MUPJM  MDPJM  KVOPJM  KVUPJM  KVDPJM  (PJM) 

MOPJ1  MUPJ1  MDPJ1   KVOPJ 1  KVUPJ 1  KVDPJ 1  (PJ1) 

MOPJ1M  MUPJ1M  MDPJ1M  KVOP1M  KVUP1M  KVDP1M  (PJ1M) 

MOPEB  MUPEB  MDPEB  KVOPEB  KVUPEB  KVDPEB  (PEB) 

MOPEBM  MUPEBM  MDPEBM  KVOPBM  KVUPBM  KVDPBM  (PEBM) 

MOPTEB  MUPTEB  MDPTEB  KVOPTB  KVUPTB  KVDPTB  (PTEB) 

MOPTBM  MUPTBM  MDPTBM  KVOBTM  KVUBTM  KVDBTM  (PTEBM) 
INPUT/ OUTPUT: 

INPUT  IS  TAKEN  FROM  A  FILE  CREATED  FROM  A  TAPE 
PROVIDED  BY  NPRDC 

OUTPUT  IS  TO  THE  PRINTER 


**      FIXED  PARAMETERS 

PARAMETER  (MXYR=7,  MXLOS=31,  MXGRD=10,  MXMOS=15) 
**       INPUT  PARAMETERS 

INTEGER  ST1,  ST2,  LYR 

INTEGER  SLOS1,  SLOS2 

INTEGER  SMOS(MXMOS),  NMOS 

INTEGER  SGRD(MXGRD),  NGRD 
**       VARIABLES  IN  ORIGINAL  SPACE 

INTEGER  T,  FLAG,  IVYR(MXYR) 

INTEGER*2  D(MXMOS,  MXLOS,  MXGRD) 

INTEGER  KV(MXYR) 

REAL  AVCINV,AINV,ATT,ATT1 

REAL  CINV(MXMOS, MXLOS, MXGRD, MXYR) ,Y(MXMOS , MXLOS, MXGRD, MXYR) 

REAL  AC INV(MXMOS, MXLOS , MXGRD) 

REAL  PJ(MXMOS, MXLOS, MXGRD),  PJM(MXMOS, MXLOS, MXGRD) 

REAL  PMLE (MXMOS, MXLOS, MXGRD),  PJLT(MXMOS , MXLOS, MXGRD) 

REAL  PJ1( MXMOS, MXLOS, MXGRD),  PJ1M(MXM0S, MXLOS, MXGRD) 

REAL  PEB (MXMOS, MXLOS, MXGRD),  PE BM( MXMOS , MXLOS, MXGRD) 

REAL  PTEB (MXMOS, MXLOS, MXGRD),  PTEBM( MXMOS , MXLOS , MXGRD) 

REAL  RRTSCA(MXYR) ,RRJ(MXYR) ,RRJLT(MXYR) .RRML(MXYR) 

REAL  RRTSAM(MXYR) ,RRJM(MXYR) ,RRJ1(MXYR) ,RRJ1M(MXYR) 

REAL  RRPEB(MXYR),RRPEBM(MXYR) 

REAL  RRPTE B ( MXYR ),RRPTBM( MXYR) 
**       VARIABLES  IN  TRANSFORMED  SPACE 

REAL  X( MXMOS, MXLOS, MXGRD, MXYR) 

REAL  XT( MXMOS, MXLOS, MXGRD, MXYR) 

REAL  XB( MXMOS, MXLOS, MXGRD),  TSCA( MXMOS, MXLOS, MXGRD) 

REAL  XTB( MXMOS, MXLOS, MXGRD) 

REAL  X J( MXMOS , MXLOS , MXGRD ) ,  X JLT( MXMOS , MXLOS , MXGRD ) 

REAL  XEB( MXMOS, MXLOS, MXGRD),  XJ 1( MXMOS , MXLOS, MXGRD) 

REAL  XTEB( MXMOS, MXLOS, MXGRD) 

REAL  XMLE( MXMOS, MXLOS, MXGRD) 

REAL  SHREB( MXMOS, MXLOS, MXGRD) 

REAL  SHRTEB( MXMOS, MXLOS, MXGRD) 

REAL  RML(MXYR),  RSL( MXYR) ,RXEB( MXYR) 

REAL  RXTEB(MXYR) 

REAL  RTSCA( MXYR), RTJ( MXYR), RTJLT( MXYR) 

REAL  R J( MXYR ),RJLT( MXYR),  RJl(MXYR) 

REAL  TSCAM( MXMOS, MXLOS, MXGRD) 
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irk 


REAL  V(MXMOS,MXLOS,MXGRD) 
REAL  VT(MXMOS,MXLOS,MXGRD) 
REAL*8  SHR,  SHRX,  SHR1,  SSB,  SST,  SSE,  XBB 
VARIABLES  TOR  FIGURE  Of  MERIT 


INTEGER  KVOTSA(MXYR) ,KVOTSM(MXYR) ,KVOPJ(MXYR) ,KVOPJl(MXYR) 

INTEGER  KVOPJM(MXYR) ,KVOPEB(MXYR) ,KV0P1M(MXYR) ,KVOPBM(MXYR) 

INTEGER  KVOPTB(MXTR) ,KVOBTM(MXYR) 

INTEGER  KVOPML(MXYR) 

INTEGER  KVUTSA(MXYR) ,KVUTSM(MXYR) ,KVUPJ(MXYR) ,KVUPJ1(MXYR) 

INTEGER  KVUPJM(MXYR) ,KVUPEB(MXYR) ,KVUP1M(MXYR) ,KVUPBM(MXYR) 

INTEGER  KVUPTB( MXYR ),KVUBTM( MXYR) 

INTEGER  KVUPML(MXYR) 

INTEGER  KVDTSA(MXYR) ,KVDTSM(MXYR) ,KVDPJ(MXYR) ,KVDPJ1(MXYR) 

INTEGER  KVDPJM(MXYR) ,KVDPEB(MXYR) ,KVDP1M(MXYR) ,KVDPBM(MXYR) 

INTEGER  KVDPTB( MXYR), KVDBTM( MXYR) 

INTEGER  KVDPML(MXYR) 

REAL  MUTSCA(MXYR) ,MUTSAM(MXYR) ,MUPJ(MXYR) ,MUPJ1(MXYR) ,MUPJM(MXYR) 

REAL  MUPJIM(MXYR) ,MUPEBM(MXYR) ,MUPEB(MXYR) ,MUPMLE(MXYR) 

REAL  MUPTBMC MXYR ) , MUPTEB ( MXYR ) 

REAL  MOTSCA(MXYR) ,MOTSAM(MXYR) ,MOPJ(MXYR) ,MOPJl(MXYR) , MOP JM( MXYR) 

REAL  MOPJIM(MXYR) ,MOPEBM(MXYR) ,MOPEB(MXYR) ,MOPMLE(MXYR) 

REAL  MOPTBM( MXYR ),MOPTEB( MXYR) 

REAL  MDTSCA(MXYR) ,MDTSAM(MXYR) ,MDPJ(MXYR) .MDPJl(MXYR) ,MDPJM(MXYR) 

REAL  MDPJIM(MXYR) ,MDPEBM(MXYR) .MDPEB(MXYR) ,MDPMLE(MXYR) 

REAL  MD  PTBM( MXYR ),MDPTEB( MXYR) 

REAL  MILLER 

**       INPUT  DATA  AREAS 

INTEGER  MOS,  LOS,  GRADE,  YR,  INV 
INTEGER  MOS 1, LOS  1, GRADE 1,YR1,INV1 

**       INPUT  INITIALIZATION 
DATA  KV/MXYR*0/ 

DATA  MOPMLE/MXYR*0/ ,MUPMLE/MXYR*0/ ,MDPMLE/MXYR*0/ 
DATA  KVOPML/MXYR*0/ ,KVUPML/MXYR*0/ ,KVDPML/MXYR*0/ 
DATA  MOTSCA/MXYR*0/ ,MUTSCA/MXYR*0/ ,MDTSCA/MXYR*0/ 
DATA  KVOTSA/MXYR*0/  ,KVUTSA/MXYR*0/  ,KVDTSA/MXYR*0/ 
DATA  MOTSAM/MXYR*0/ ,MUTSAM/MXYR*0/ ,MDTSAM/MXYR*0/ 
DATA  KVOTSM/MXYR*0/ ,KVUTSM/MXYR*0/ ,KVDTSM/MXYR*0/ 
DATA  MOPJ/MXYR*0/ ,MUPJ/MXYR*0/ ,MDPJ/MXYR*0/ 
DATA  KVOPJ/MXYR*0/ ,KVUPJ/MXYR*0/ ,KVDPJ/MXYR*0/ 
DATA  MOPJM/MXYR*0/  ,MUPJM/MXYR*0/  ,MDPJM/MXYR*0/ 
DATA  KVOPJM/MXYR*0/ ,KVUPJM/MXYR*0/ ,KVDPJM/MXYR*0/ 
DATA  MOPJl/MXYR*0/ ,MUPJ1/MXYR*0/ ,MDPJ1/MXYR*0/ 
DATA  KVOPJl/MXYR*0/ ,KVUPJ1/MXYR*0/ ,KVDPJ1/MXYR*0/ 
DATA  MOPJ1M/MXYR*0/  ,MUPJ1M/MXYR*0/  ,MDPJ1M/MXYR*0/ 
DATA  KV0P1M/MXYR*0/ ,KVUP1M/MXYR*0/ ,KVDP1M/MXYR*0/ 
DATA  MOPEB/MXYR*0/ ,MUPEB/MXYR*0/ ,MDPEB/MXYR*0/ 
DATA  KVOPEB/MXYR*0/ ,KVUPEB/MXYR*0/ ,KVDPEB/MXYR*0/ 
DATA  MOPEBM/MXYR*0/ ,MUPEBM/MXYR*0/ ,MDPEBM/MXYR*0/ 
DATA  KVOPBM/MXYR*0/ , KVUPBM/MXYR*0/ , KVDPBM/MXYR*0/ 
DATA  MOPTEB/MXYR*0/ , MUPTEB/ MXYR* 0/ ,MDPTEB/MXYR*0/ 
DATA  KVOPTB/MXYR*0/,KVUPTB/MXYR*0/,KVDPTB/MXYR*0/ 
DATA  MOPTBM/MXYR*0/ ,MUPTBM/MXYR*0/ ,MDPTBM/MXYR*0/ 
DATA  KVOBTM/MXYR*0/ ,KVUBTM/MXYR*0/ ,KVDBTM/MXYR*0/ 
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**       READ  PARAMETER  DATA  TO  USE  IN  THIS  RUN.  ** 


100 


101 


102 


103 


104 


105 


WRITE(5.100) 

FORMAT(f  ENTER  1ST  YEAR  AND  LAST  YEAR  TO  USE  FOR  ESTIMATION') 

READ(5,*)  ST1,ST2 

WRITE(6,*)  'ESTIMATION  YEARS:  ',ST1,ST2 

WRITE(5,101) 

FORMAT(f  ENTER  LAST  YEAR  PRESENT  IN  THE  DATA  BASE*) 

READ(5,*)  LYR 

WRITE(6,*)  'LAST  YEAR  AVAILABLE  FOR  VALIDATION:  '  ,LYR 

WRITE(5.102) 

FORMAT(f  ENTER  NO.  OF  MOS  FOLLOWED  BY  ARRAY  OF  SAME  LENGTH') 

READ(5,*)  NMOS,  (SMOS(I),  1=1, NMOS) 

WRITE(6,*)  'MOS  SELECTED:',  (SMOS(I),  1=1, NMOS) 

WRITE(5.103) 

FORMAT(f  ENTER  1ST  AND  LAST  LOS  VALUE  TO  USE') 

READ(5,*)  SLOS1.SLOS2 

WRITE(6,*)  'LOS  RANGE:  ' ,SLOSl,SLOS2 

WRITE(5.104) 

FORMAT(f  ENTER  NO.  OF  GRADE  FOLLOWED  BY  ARRAY  OF  SAME  LENGTH') 

READ(5,*)  NGRD,  (SGRD(I),  1=1, NGRD) 

WRITE(6,*)  'GRADES  SELECTED',  (SGRD(I),  1=1, NGRD) 

WRITE(5.105) 

FORMAT(f  ENTER  DEE  FACTOR') 

READ(5,*)  DEE 

WRITE(6,*)  'DEE  FACTOR: ', DEE 


COMPUTE  START  AND  END  YEARS  FOR  VALIDATION 
NYR=LYR-ST1+1 
NVYR=LYR-ST2 
NLOS=SLOS2-SLOSl+l 
T=ST2-ST1+1 
DO  41  I=1,NVYR 
IVYR(I)=ST2+I 
41    CONTINUE 

**       INITIALIZE  D  ARRAY  TO  ALLOW  FINDING  MISSING  VALUES.        ** 

DO  1  I=l,MXMOS 
DO  2  J=l,MXLOS 
DO  3  K=1,MXGRD 
D(I,J,K)=-9999 
3       CONTINUE 
2      CONTINUE 
1     CONTINUE 

**      TAPE  PROCESSING  TO  INPUT  INVENTORY  AND  ATTRITION  VALUES,   ** 
**       DEVELOPE  D  MATRIX  ACCORDING  TO  EXTRACTION  CRITERIA.        ** 

FLAG=0 
IEOF=0 
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* 

11 


** 


12 


8 


CALL  READER(M0S1, LOS  1, GRADE 1 ,YR1 , INV1,ATT1,  NMOS,SMOS, 

*  NGRD,SGRD,  SL0S1,SL0S2,  ST1.ST2,  IM1 ,IL1,IG1,IT1,IE0F) 

IF(IEOF  .NE.  0)  THEN 

WRITE(6,*)  *****  NO  DATA  MEETS  SELECTIONS  REQS' 

STOP 
END  IF 

IF(IEOF  .NE.  0)  GO  TO  9 
MOS=MOSl 
LOS=LOSl 
GRADE=GRADE1 
YR=YR1 
INV=INV1 
ATT=ATT1 
IM=IM1 
IL=IL1 
IG=IG1 
IT=IT1 
IF(YR.  GT.  ST1)  FLAG  =  1 

CHECK  CASE  WHERE  DATA  BEGINS  PAST  1ST  YR 
IF(YR  .GT.  ST1)  CINV(IM,IL,IG,IT-1)=. 5*FLOAT( INV) 
CALL  READER(M0S1,L0S1,GRADE1,YR1,INV1,ATT1,  NMOS,SMOS, 

*  NGRD,SGRD,  SLOSl,SLOS2,  ST1,ST2,  IM1,IL1,IG1,IT1,IE0F) 
ONLY  FOR  CHECKING  PURPOSES.  MARK  D  TO  INDICATE  SOME  YR  PRESENT 

D(IM,IL,IG)=0 

IF(.NOT.  (IEOF.  EQ.  0  .AND.  MOS1.EQ.  MOS  .AND.  LOS1.  EQ.  LOS 

*  .AND.  GRADE1.EQ.  GRADE))  GO  TO  8 
CENTRAL  INV.  FOR  YR1  -  1 

IF(YR1  .GT.  YR+1)  THEN 

CINV( IM1 , IL1 , IG1 , IT1-1)=. 5*FLOAT( INV1) 

AINV=. 5*FL0AT( INV) 

I=YR+1 

FLAG=1 
ELSE 

AINV=. 5*FLOAT( INV+INV1) 
END  IF 

CINV(IM,IL,IG,IT)=AMAX1(AINV,ATT) 
Y(IM,IL,IG,IT)=ATT 
YR=YR1 
INV=INV1 
ATT=ATT1 
IM=IM1 
IL=IL1 
IG=IG1 
IT=IT1 
CALL  READER(M0S1,L0S1,GRADE1,YR1,INV1,ATT1,  NMOS,SMOS, 

*  NGRD.SGRD,  SLOSl,SLOS2,  ST1.ST2,  IM1,IL1 , IG1.IT1, IEOF) 
GO  TO  12 

CONTINUE 

WHEN  YEARS  MISSING  AT  THE  END 
IF(YR  .LT.  LYR)  THEN 

FLAG=1   - 

AINV=.  5*FLOAT( INV) 
ELSE 
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AINV=FLOAT(INV) 
END  IF 

CINV(IM,IL,IG,IT)=AMAX1(AINV,ATT) 
Y(IM,IL,IG,IT)=ATT 

CC=CINV(IM,IL,IG,IT) 
TEMP=-1.  +2.  *ATT/(  1.  +CC) 
TEMP1=-1.  +  2.*(1.+ATT)/(1. +CC) 
CC=ABS(TEMP) 
CC1=ABS( TEMPI) 

IFCCC   .  GT.    1.    .OR.    CC1   .  GT.    1.  )  VRITE(6,*)    '***  TEMP,   TEMP1=' , 
*  TEMP, TEMPI,' ***MOS, LOS, GR,YR=' ,SMOS( IM) ,IL,SGRD( IG) ,IT 

GO  TO   11 

CONTINUE 

k  A  A  AAA********irA*-A-A-A-ArA*-A*-A*-A-A*-A-Ai 


i   t   ii  r      1 .  «   t   t  t      *      t   i   t   t   t   i   i   t   i   i   i__i._i   i..i.  i   t   i   i 

nv^nnnnnnnr  TC7c"~"c7r"io?7C7r7r7f  7"7vTr"7orTr7r 

CALCULATE  AVERAGE  CENTRAL  INVENTORY  AND  MAXIMUM  ** 

LIKELIHOOD  ESTIMATE.  ** 

DO  10  IM=l,NMOS 

DO  20  IL=l,NLOS 
DO  30  IG=1,NGRD 
SUMY=0 

ACINV(IM,IL,IG)=0 
DO  35  IT=1,T 
SUMY=SUMY+Y( IM , IL , IG , IT) 

ACINV( IM, IL, IG)=ACINV( IM, IL, IG)+CINV( IM, IL, IG, IT) 
CONTINUE 
IF(ACINV(IM,IL,IG)  .  NE.  0)  THEN 

PMLE( IM, IL, IG)=SUMY/ACINV( IM, IL, IG) 
ELSE 

PMLE(IM,IL,IG)=0 
ENDIF 

ACINV(IM,IL,IG)=ACINV(IM,IL,IG)/T 
DO  40  IT=1,NYR 
IF  ANY  CENTRAL  INVENTORY>0  (OVER  EST.  YRS  T)  THEN  D=l 

IF(CINV(IM,IL,IG,IT).GT.  0  .AND.  IT.  LE.  T)  D(  IM,  IL,  IG)=1 
IF  ANY  CENTRAL  INVENTORY>0  (OVER  VAL.  YRS  >  T)  THEN  COUNT 

IF(CINV(IM,IL,IG,IT).GT. 0  .AND.  IT. GT. T)  KV( IT-T)=KV(IT-T)+1 
TRANSFORMATION  OF  CINV  USING  THE  FREEMAN-TUKEY 
DOUBLE  ARCS IN  FORMULA. 
CALL  FTT(CINV(IM,IL,IG,IT),Y(IM,IL,IG,IT),X(IM,IL,IG,IT)) 
XT(IM,IL,IG,IT)=X(IM,IL,IG,IT)   /   SQRT(0.  5+CINV( IM, IL,IG, IT) ) 
CONTINUE 
CONTINUE 
CONTINUE 
CONTINUE 


r-AA-AAiti? 


:1rtfic!;i?irfrifAiTititiricAifAic. 


'i-.Wr  vWf  irAir  A^cAA  A  -A  A A  A A  A  ■A,;'t  A  A  A  A  A  A  AAirAirirAA 

CHECK  FOR  MISSING  DATA,  COMPUTE  THE  TRANSFORMED  ** 

CELL  AVERAGES  AND  GRAND  MEAN.  ** 

^AAAAAAAAAAAAAAAAAAAirAir-AirAirAirMcirAir>tAAirAA-Ai(AAAAAAAAA-/iAAAAA*AA-!tA 

XBB=0 
KK=0 

DO  50   IM=l,NMOS  ** 
DO  60   IL=l,NLOS 
DO  70  IG=1,NGRD 
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*rt 


IF(D(IM,IL,IG)  .EQ.  -9999)  THEN 

REPORT  MISSING  COMBINATION  IM,IL,IG  &  CLEAR  D  ENTRY 
D(IM,ILsIG)=0 
ILD=SLOSl+IL-l 
FLAG=1 
END  IF 

KK=KK+D(IM,IL,IG) 
XB(IM,IL,IG)=0 
XTB(IM,IL,IG)=0 
DO  80  IT=1,T 
XB(IM,IL,IG)=XB(IM,IL,IG)+X(IM,IL,IG,IT) 
XTB( IM , IL, IG)=XTB( IM, IL, IG)+XT( IM , IL, IG , IT) 
CONTINUE 

XB(IM,IL,IG)=XB(IM,IL,IG)/T 
XTB( IM, IL, IG)=XTB( IM , IL, IG)/T 
EMPIRICAL  BAYES  PREPARATION  ON  THE  TRANSFORMED  SCALE. 
XEB(IM,IL,IG)=XB(IM,IL,IG) 
XTEB(IM,IL,IG)=XTB(IM,IL,IG) 

XBB=XBB+D( IM, IL, IG)*XB( IM, IL, IG) 

CONTINUE 
CONTINUE 
CONTINUE 
XBB=XBB/KK 

**       BAYES  VARIANCE  ITERATION  SECTION.  ** 

IC=0 

COMPUTE  TRANSFORMED  SCALE  CELL  VARIANCE. 
DO  500  IM=1,NM0S 


80 


** 


70 
60 
50 


499 


522 


520 

510 

500 
* 

530 


533 
532 
531 


CINV(IM,IL,IG,IT)) 
*  1.5708 


DO  510  IL=1,NL0S 
DO  520  IG=1,NGRD 

IF  (D(IM,IL,IG).NE.  1)  GOTO  520 
AVSQR  =  0.  0 
DO  522  IT=1,T 

AVSQR  =  AVSQR  +  SQRT(0.5  + 
Z  =  XEB(IM,IL,IG)  +  (AVSQR/T) 
V(IM,IL,IG)  =  VAR(Z) 
CONTINUE 
CONTINUE 
CONTINUE 


A=0 
AI=A 
SALPH=0 

DO  531  IM=1,NM0S 
DO  532  IL=1,NL0S 
DO  533  IG=1,NGRD 

IF  (D(IM,IL,IG).EQ.  1)  SALPH  =  SALPH  +  1. /(A  +  V(IM,IL,IG)) 
CONTINUE 
CONTINUE 

ZB=0 

COMPUTE  BAYES'  ALPHA  AND  GAMMA  PARAMETERS 
DO  537  IM=1,NM0S 
DO  538  IL=1,NL0S 


MAR03850 
MAR03860 
MAR03870 
MAR03880 
MAR03890 
MAR03900 
MAR03910 
MAR03920 
MAR03930 
MAR03940 
MAR03950 
MAR03960 
MAR03970 
MAR03980 
MAR03990 
MAR04000 
MAR04010 
MAR04020 
MAR04030 
MAR04040 
MAR04050 
MAR04060 
MAR04070 
MAR04080 
MAR04090 
MAR04100 
MAR04110 
MAR04120 
MAR04130 
MAR04140 
MAR04150 
MAR04160 
MAR04170 
MAR04180 
MAR04190 
MAR04200 
MAR04210 
MAR04220 
MAR04230 
MAR04240 
MAR04250 
MAR04260 
MAR04270 
MAR04280 
MAR04290 
MAR04300 
MAR04310 
MAR04320 
MAR04330 
MAR04340 
MAR04350 
MAR04360 
MAR04370 
MAR04380 
MAR04390 
MAR04400 
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DO  539  IG=1,NGRD 

IF  (D(IM,IL,IG).EQ.  1)  THEN 
ALPH  =  l./(A  +  V(IM,IL,IG)) 
GAM  =  ALPH  /  SALPH 
ZB   =  ZB  +  GAM*XB(IM,IL,IG) 
END  IF 
CONTINUE 
CONTINUE 
CONTINUE 
F=0 
G=0 
DO  540  IM=1,NM0S 
DO  541  IL=1,NL0S 
DO  542  IG=1,NGRD 

IF  (D(IM,IL,IG).EQ. 1)  THEN 


/(A  +  V(IM,IL,IG)) 

+  (ALPH  *  (XB(IM,IL,IG)  -  ZB)**2) 

+  (ALPH  *  (XB(IM,IL,IG)-ZB))**2 


ALPH  =  1. 
F    =  F 
G    =  G 
ENDIF 
CONTINUE 
CONTINUE 
CONTINUE 

COMPUTE  THE  PRIOR  VARIANCE  (A) 
A  =  A  -  (KK-1-F)/G 

IF  (A.  LE.  0.  )  THEN 
A=0 

GOTO  550 
ENDIF 
IF  (ABS(A-AI).GT.  .0001)  GOTO  530 

COMPUTE  THE  EMPIRICAL  BAYES  ATTRITION  RATES 
DO  570  IM=1,NM0S 
DO  580  IL=1,NL0S 
DO  590  IG=1,NGRD 

IF  (D(IM,IL,IG).EQ.  1)  THEN 
DEN  =  A  +  V(IM,IL,IG) 

XEB(IM,IL,IG)  =  (A/DEN)*XB(IM,IL,IG)+((DEN-A)/DEN  )*ZB 
ENDIF 
CONTINUE 
CONTINUE 
CONTINUE 
IC=IC+1 

IF  (A.  EQ.  0)  GOTO  599 
IF  (IC.LT.  10)  GOTO  499 

CONTINUE 

COMPUTE  THE  EMPIRICAL  BAYES  SHRINKAGE  RATES 
DO  600  IM=1,NM0S 
DO  601  IL=1,NL0S 
DO  602  IG=1,NGRD 

SHREB(IM,IL,IG)  =  0 

IF  (D(IM,IL,IG).EQ.  1)  THEN 

SHREB(IM,TL,IG)  =  V( IM,IL,IG)/(A  +  V(IM,IL,IG)) 
ENDIF 
CONTINUE 


MAR04410 
MAR04420 
MAR04430 
MAR04440 
MAR04450 
MAR04460 
MAR04470 
MAR04480 
MAR04490 
MAR04500 
MAR04510 
MAR04520 
MAR04530 
MAR04540 
MAR04550 
MAR04560 
MAR04570 
MAR04580 
MAR04590 
MAR04600 
MAR04610 
MAR04620 
MAR04630 
MAR04640 
MAR04650 
MAR04660 
MAR04670 
MAR04680 
MAR04690 
MAR04700 
MAR04710 
MAR04720 
MAR04730 
MAR04740 
MAR04750 
MAR04760 
MAR04770 
MAR04780 
MAR04790 
MAR04800 
MAR04810 
MAR04820 
MAR04830 
MAR04840 
MAR04850 
MAR04860 
MAR04870 
MAR04880 
MAR04890 
MAR04900 
MAR04910 
MAR04920 
MAR04930 
MAR04940 
MAR04950 
MAR04960 
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601     CONTINUE 
600   CONTINUE 

**       ALTERNATE  BAYES  VARIANCE  ITERATION  SECTION.  ** 


699 


725 


720 
710 
700 

730 


733 
732 
731 


739 
738 
737 


COMPUTE  TRANSFORMED  SCALE  CELL  VARIANCE. 
DO  700  IM=1,NM0S 
DO  710  IL=1,NL0S 
DO  720  IG=1,NGRD 
VT(IM,IL,IG)  =0.0 
IF  (D(IM,IL,IG).NE. 1)  GOTO  720 
DO  725  IT=1,T 

ZT  =  X(IM,IL,IG,IT)  +  1.5708  *  SQRT(0.5  + 

*  CINV(IM,IL,IG,IT)) 

VT(IM,IL,IG)  =  VT(IM,IL,IG)  +  VAR(ZT)/(0.5  + 

*  CINV(IM,IL,IG,IT)) 
VT(IM,IL,IG)  =  VT(IM,IL,IG)/T**2 

CONTINUE 
CONTINUE 
CONTINUE 

AT=0 
AT I = AT 
SALPHT=0 

DO  731  IM=1,NM0S 
DO  732  IL=l,NLOS 
DO  733  IG=1,NGRD 

IF  (D(IM,IL,IG).EQ.  1)  SALPHT=  SALPHT+  1. /(AT+  VT(IM,IL,IG)) 
CONTINUE 
CONTINUE 
ZTB=0 

COMPUTE  BAYES  ALPHA  AND  GAMMA  PARAMETERS 
DO  737  IM=l,NMOS 
DO  738  IL=l,NLOS 
DO  739  IG=1,NGRD 

IF  (D(IM,IL,IG).EQ.  1)  THEN 
ALPHT=  1. /(AT+  VT(IM,IL,IG)) 
GAMT  =  ALPHT  /  SALPHT 
ZTB  =  ZTB+  GAMT*XTB(IM,IL,IG) 
ENDIF 
CONTINUE 
CONTINUE 
CONTINUE 

PZTB=0.  5*(  1. +SIN(  ZTB)  ) 
FT=0 
GT=0 

DO  740  IM=l,NMOS 
DO  741  IL=l,NLOS 
DO  742  IG=1,NGRD 

IF  (D(IM,IL,IG).EQ.  1)  THEN 
ALPHT= 
FT   = 
GT 


1. /(AT  +  VT(IM,IL,IG)) 

FT  +  (ALPHT  *  (XTB(IM,IL,IG) 


ZTB)**2) 
=  GT  +  (ALPHT  *  (XTB(IM,IL,IG)  -  ZTB))**2 


742 
741 


ENDIF 
CONTINUE 
CONTINUE 


MAR04970 
MAR04980 
MAR04990 
MAR05000 
MAR05010 
MAR05020 
MAR05030 
MAR05040 
MAR05050 
MAR05060 
MAR05070 
MAR05080 
MAR05090 
MAR05100 
MAR05110 
MAR05120 
MAR05130 
MAR05140 
MAR05150 
MAR05160 
MAR05170 
MAR05180 
MAR05190 
MAR05200 
MAR05210 
MAR05220 
MAR05230 
MAR05240 
MAR05250 
MAR05260 
MAR05270 
MAR05280 
MAR05290 
MAR05300 
MAR05310 
MAR05320 
MAR05330 
MAR05340 
MAR05350 
MAR05360 
MAR05370 
MAR05380 
MAR05390 
MAR05400 
MAR05410 
MAR05420 
MAR05430 
MAR05440 
MAR05450 
MAR05460 
MAR05470 
MAR05480 
MAR05490 
MAR05500 
MAR05510 
MAR05520 


60 


740 

•kit 


750 


790 
780 

770 


802 
801 
800 


CONTINUE 

COMPUTE  THE  PRIOR  VARIANCE  (A) 
AT=  AT-  (KK-1-FT)/GT 

IF  (AT.LE.  0.  )  THEN 
AT=0 

GOTO  750 
ENDIF 
IF  (ABS(AT-ATI).GT.  .0001)  GOTO  730 

COMPUTE  THE  EMPIRICAL  BAYES  ATTRITION  RATES 
DO  770  IM=1,NM0S 
DO  780  IL=l,NLOS 
DO  790  IG=1,NGRD 

IF  (D(IM,IL,IG).EQ.  1)  THEN 
DENT  =  AT  +  VT(IM,IL,IG) 

XTEB(IM,IL,IG)  =  (AT/DENT)*XTB(IM,IL,IG)  + 
"  ((DENT-AT)/DENT)*ZTB 

ENDIF 
CONTINUE 
CONTINUE 
CONTINUE 

COMPUTE  THE  EMPIRICAL  BAYES  SHRINKAGE  RATES 
DO  800  IM=l,NMOS 
DO  801  IL=l,NLOS 
DO  802  IG=1,NGRD 

SHRTEB(IM,IL,IG)  =  0 

IF  (D(IM,IL,IG).EQ.  1)  THEN 

SHRTEB(IM,IL,IG)  =  VT( IM,IL,IG)/(AT  +  VT(IM,IL, IG)) 
ENDIF 
CONTINUE 
CONTINUE 
CONTINUE 


** 


**       JAMES -STEIN  SHRINKAGE  RATES.  ** 

SST=0 
SSE=0 

SSB=0 

COMPUTE  THE  SUMS  OF  SQUARES  (SST,  SSE  AND  SSB) 
DO  90  IM=1,NM0S 
DO.  95  IL=1,NL0S 
DO  110  IG=1,NGRD 
TEMP 1=0 
TEMP2=0 
DO  120  IT=1,T 
TEMP1=TEMP1+(X( IM, IL, IG, IT) -XBB)**2 
TEMP2=TEMP2+(X( IM, IL, IG, IT) -XB( IM, IL, IG) )**2 
CONTINUE 

SST=SST+D( IM , IL , IG)*TEMP1 
SSE=SSE+D(IM,IL,IG)*TEMP2 
CONTINUE 
CONTINUE 
CONTINUE 
SSB=SST-SSE   - 

COMPUTE  THE  JAMES-STEIN  SHRINKAGE  RATES 
SHRX=(SSE/SSB)*  (KK-3)/(KK*(T-l)+2) 


120 


110 

95 

90 

** 


MAR05530 
MAR05540 
MAR05550 
MAR05560 
MAR05570 
MAR05580 
MAR05590 
MAR05600 
MAR05610 
MAR05620 
MAR05630 
MAR05640 
MAR05650 
MAR05660 
MAR05670 
MAR05680 
MAR05690 
MAR05700 
MAR05  710 
MAR05720 
MAR05730 
MAR05  740 
MAR05  750 
MAR05760 
MAR05770 
MAR05780 
MAR05790 
MAR05800 
MAR05810 
MAR05820 
MAR05830 
MAR05840 
MAR05850 
MAR05860 
MAR05870 
MAR05880 
MAR05890 
MAR05900 
MAR05910 
MAR05920 
MAR05930 
MAR05940 
MAR05950 
MAR05960 
MAR059  70 
MAR05980 
MAR05990 
MAR06000 
MAR06010 
MAR06020 
MAR06030 
MAR06040 
MAR06050 
MAR06060 
MAR06070 
MAR06080 
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SHR=DMIN1(SHRX,  1.  DO) 
SHR=DMAX1(SHR  ,  O.DO) 
SHR1=T*SHR 

t.-r-.i    _t      i.    t       i    _t. _!__*-    i..i.    11       t       i       t    ,t      i      1,.1,-t..!-,  I .,  t  ulii  I  ii  |„',,1 ,,  1..  1..!..f      tittitti    _t_.t..t_.t    .1.    i.    i..i      i 

**  LIMITED  TRANSLATION  OPTION  &  JAMES-STEIN 

**  ON  TRANSFORMED  SCALE. 

NN=KK*(T-1) 
DO   125    IM=l,NMOS 
DO   130   IL=1,NL0S 
DO   140  IG=1,NGRD 
XJ(IM,IL,IG)=0 
XJ1(IM,IL,IG)=0 
IF  (D(IM,IL,IG).EQ.  1)  THEN 

XJ(IM,IL,IG)=XBB+(1-SHR)*(XB(IM,IL,IG)-XBB) 

XJ1(IM,IL,IG)=XBB+(1-SHR1)*(XB(IM,IL,IG)-XBB) 

XJLT(IM,IL,IG)=XBB+(1=SHR*RH0(XB(IM,IL,IG),XBB,SSE, 

SSB,DEE,KK,T,NN))*(XB(IM,IL,IG)-XBB) 
XMLE( IM, IL, IG)=SQRT( .  5+ACINV( IM, IL , IG) ) 
*  *ASIN(2*PMLE(IM,IL,IG)-1) 

ENDIF 
CONTINUE 
CONTINUE 
CONTINUE 


140 
130 
125 

»f,-*-  i   1   I   t   I  _t   '■.l.,l„'.J„L.1.,l,.t-t..1..l..l..l   «  -t.  I..1   t   I   I   |   t  .t.i.t  .1   t   ♦__>._» — I 1 | I  ..  I  .  .1  .  .  t  ..  I  .t„l„»--r„1„t,-f-,L.',,JM«l—T«.T«  .!,.1M,|L»'«  .'-.!■  ■'■■la»ta»la 


**  VALIDATION.      RISK  EVALUATED   IN  TRANSFORM  SPACE. 

AA  AA1  Ari^r^r^*:^yrfrVrycyr*?Wrftyc^ 


** 


c  *** 


C  *** 


180 
170 
160 


DO   150   IT=1,   NVYR 
RTSCA( IT)=0 
RJ(IT)=0 
RJ1(IT)=0 
RXEB(IT)=0 
RXTEB(IT)=0 
RJLT(IT)=0 
RML(IT)=0 
DO   160   IM=1,NM0S 
DO   170   IL=1,NL0S 
DO   180   IG=1,NGRD 
IF(D(IM,IL,IG).EQ.  1   .AND.    CINV(  IM,IL,IG,IT+T).  NE.  0)   THEN 
RTSCA( IT)=RTSCA( IT)+(X( IM, IL, IG ,T+IT) -     XB( IM, IL, IG) )**2 
RJ(IT)=       RJ(IT)+(X(IM,IL,IG,T+IT)-     XJ(IM,IL,IG) )**2 
RJ1(IT)=  RJ1(IT)+(X(IM,IL,IG,T+IT)-   XJ1(IM, IL,IG))**2 
RXEB( IT)=RXEB( IT)+(X( IM, IL, IG,T+IT) -   XEB( IM, IL, IG) )**2 
RXTEB(IT)=RXTEB(IT)+((X(IM,IL,IG,T+IT)/SQRT(0.5  + 

CINV(IM,IL,IG,IT+T))-   XTEB( IM,IL,IG))**2) 
/VT(IM,IL,IG) 
RJLT  (IT)=  RJLT(IT)+(X(IM,IL,IG,T+IT)-XJLT(IM,IL,IG))**2 
RML     (IT)=  RML  (IT)+(X( IM,IL,IG,T+IT) -XMLE( IM,IL,IG))**2 
ENDIF 
CONTINUE 
CONTINUE 
CONTINUE 

RTSCA( IT)=RTSCA( IT) /KV( IT) 
RJ     (IT)=RJ     (IT)/KV(IT) 
RJ1   (IT)=RJ1  (IT)/KV(IT) 
RXEB( IT)=RXEB( IT) /KV( IT) 


MAR06090 
MAR06100 
MAR06110 
MAR06120 
'MAR06130 
MAR06140 
MAR06150 
MAR06160 
MAR06170 
MAR06180 
MAR06190 
MAR06200 
MAR06210 
MAR06220 
MAR06230 
MAR06240 
MAR06250 
MAR06260 
MAR06270 
MAR06280 
MAR06290 
MAR06300 
MAR06310 
MAR06320 
MAR06330 
MAR06340 
MAR06350 
MAR06360 
MAR06370 
MAR06380 
MAR06390 
MAR06400 
MAR06410 
MAR06420 
MAR06430 
MAR06440 
MAR06450 
MAR06460 
MAR06470 
MAR06480 
MAR06490 
MAR06500 
MAR06510 
MAR06520 
MAR06530 
MAR06540 
MAR06550 
MAR06560 
MAR06570 
MAR06580 
MAR06590 
MAR06600 
MAR06610 
MAR06620 
MAR06630 
MAR06640 
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RXTEB( IT)=RXTEB( IT)/(KV( IT)*T) 
C  ***  RJLT(IT)=RJLT(IT)/KV(IT) 
RML  (IT)=RML  (IT)/KV(IT) 
RTJ( IT)=RTSCA( IT)/RJ( IT) 
RTJLT( IT)=RTSCA( IT) /RJLT( IT) 

RSL( IT)=( RJ( IT) -RTSCA( IT) ) /( RJLT( IT) -RTSCA( IT) ) 
CONTINUE 


C  *** 

150 


XB(IM,IL,IG),  AVCINV) 
XJ(IM,IL,IG),  AVCINV) 


**       DO  RISK  ANALYSIS  ON  THE  ORIGINAL  SCALE.  ** 

DO  200  IM=1,NM0S 
DO  210  IL=1,NL0S 
DO  220  IG=1,NGRD 
**      AVERAGE  CENTRAL  INVENTORY  OVER  ESTIMATION  SET 
CINVS=0 
DO  230  IT=1,T 

AVCINV=AVCINV  +  CINV(IM, IL,IG,IT) 
230      CONTINUE 

AVCINV=AVCINV/T 

TSCA(IM,IL,IG)=D(IM,IL,IG)*SCINV( 
PJ(IM,IL,IG)=D(IM,IL,IG)*SCINV( 
PJ1(IM,IL,IG)=D(IM,IL,IG)*SCINV(  XJ1( IM,IL,IG) ,  AVCINV) 
PEB(IM,IL,IG)=D(IM,IL,IG)*SCINV(  XEB( IM,IL,IG) ,  AVCINV) 
PTEB(IM,IL,IG)=D(IM,IL,IG)*0. 5*(1.  +  SIN(XTEB( IM, IL, IG))) 
C  ***    PJLT(IM,IL,IG)=D(IM,IL,IG)*SCINV(XJLT(IM,IL,IG),  AVCINV) 
IF  (D(IM,IL,IG).EQ.  1)  THEN 

TSCAM(IM,IL,IG)=MILLER(   XB( IM,IL,IG) ,  AVCINV  ) 

PJM(IM,IL,IG)=MILLER(   XJ( IM,IL,IG) ,  AVCINV  ) 

PJ1M(IM,IL,IG)=MILLER(  XJ1(IM, IL,IG) ,  AVCINV  ) 

PEBM(IM,IL,IG)=MILLER(  XEB(IM,IL,IG) ,  AVCINV  ) 

XMIL  =  XTEB(IM,IL,IG)  *  SQRT(0.5  +  AVCINV) 
PTEBM(IM,IL,IG)=MILLER(  XMIL  ,  AVCINV  ) 
END  IF 
220     CONTINUE 
210    CONTINUE 
200   CONTINUE 

**       VALIDATION.   RISK  EVALUATED  IN  ORIGINAL  SPACE.  ** 

DO  240  IT=1,  NVYR 
RRTSCA(IT)=0 
RRTSAM(IT)=0 
RRJ(IT)=0 
RRJM(IT)=0 
RRPEB(IT)=0 
RRPTEB(IT)=0 
RRPEBM(IT)=0 
RRPTBM(IT)=0 
RRJ1(IT)=0 
RRJ1M(IT)=0 
C  ***  RRJLT(IT)=0 
RRML(IT)=0 
KP=0 

DO  250  IM=1,NM0S 
DO  260  IL=l,NLOS 
DO  270  IG=1,NGRD 


MAR06650 
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'MAR06690 
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MAR06710 
MAR06720 
MAR06730 
MAR06740 
MAR06750 
MAR06760 
MAR06770 
MAR06780 
MAR06790 
MAR06800 
MAR06810 
MAR06820 
MAR06830 
MAR06840 
MAR06850 
MAR06860 
MAR06870 
MAR06880 
MAR06890 
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MAR06940 
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MAR06990 
MAR07000 
MAR07010 
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MAR07030 
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MAR07050 
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MAR07140 
MAR07150 
MAR07160 
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MAR07180 
MAR07190 
MAR07200 
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* 

* 
* 
* 
* 
* 

* 

c   * 


CC=CINV( IM , IL, IG ,T+IT) 
IF(CC  .NE. 0)  THEN 

AR=Y(IM,IL,IG,T+IT)/CC 
ELSE 

AR=0 
END  IF 

IF(D(IM,IL,IG)  .NE.  0  .AND.  PMLE( IM,IL,IG). NE. 0 
•  AND.  PMLE(IM,IL,IG).NE.  1)  THEN 
KP=KP+1 
RRML( IT)=RRML( IT)+(CC*( AR-PMLE( IM, IL, IG) )**2) 

/(PMLE(IM,IL,IG)*(1.  -  PMLE(IM,IL,IG))) 
ENDIF 
PRODUCE  CHI -SQUARE  MEASURE  OF  EFFECTIVENESS 
RRTSCA( IT)=  RRTSCA( IT)+CHISQR(D( IM, IL, IG) ,TSCA( IM, IL, IG) , 

CC,AR) 
RRTSAM( IT)=  RRTSAM( IT)+CHISQR(D( IM, IL, IG) ,TSCAM( IM, IL, IG) , 

CC.AR) 
RRJ   (IT)=  RRJ   (IT)+CHISQR(D(IMSIL,IG),PJ  (IM,IL,IG), 

CC,AR) 
RRJM  (IT)=  RRJM  (IT)+CHISQR(D( IM, IL,IG) ,PJM  (IM,IL,IG) , 

CC,AR) 
RRJ1  (IT)=  RRJ1  (IT)+CHISQR(D(IM,IL,IG),PJ1  (IM,IL,IG), 

CC,AR) 
RRJ1M  (IT)=  RRJ1M  (IT)+CHISQR(D( IM,IL,IG) ,PJ1M(IM,IL,IG) , 

CC,AR) 
RRPEB  (IT)=  RRPEB  ( IT)+CHISQR(D( IM,IL,IG) ,PEB  (IM,IL,IG), 

CC,AR) 
RRPTEB  (IT)=  RRPTEB  (IT)+CHISQR(D( IM,IL,IG) ,PTEB  (IM,IL,IG), 

CC,AR) 
RRPEBM( IT)=  RRPEBM( IT)+CHISQR(D( IM, IL, IG) ,PEBM( IM, IL, IG) , 

CC,AR) 
RRPTBM( IT)=  RRPTBM( IT)+CHISQR(D( IM, IL, IG) ,PTEBM( IM, IL, IG) , 

CC,AR) 
RRJLT  (IT)=  RRJLT  ( IT)+CHISQR(D( IM, IL,IG) ,PJLT(IM,IL,IG) , 

CC,AR) 
PRODUCE  MEAN  ABSOLUTE  DEVIATION  MEASURE  OF  EFFECTIVENESS 
IF(D(IM,IL,IG)  .NE.  0  .AND.  PMLE(  IM,IL,IG).  NE.  0 
.AND.  PMLE(IM,IL,IG).NE.  1)  THEN 
CALL  UOM(  CC,  AR,  PMLE( IM,IL,IG) ,  MOPMLE(IT), 

MUPMLE(IT),  KVOPML(IT),  KVUPML( IT) ,  KVDPML(IT)) 
ENDIF 

IF(CC.NE.  0.  )  THEN 

CALL  UOM(  CC,  AR,  TSCA(IM,IL,IG) ,  MOTSCA(IT), 

MUTSCA(IT),  KVOTSA(IT),  KVUTSA( IT) ,  KVDTSA(IT)) 

CALL  UOM(  CC,  AR,  TSCAM(IM, IL.IG) ,  MOTSAM(IT), 

MUTSAM(IT),  KVOTSM(IT),  KVUTSM( IT) ,  KVDTSM(IT)) 

CALL  UOM(  CC,  AR,  PJ( IM,IL,IG) ,  MOPJ(IT), 

MUPJ(IT),  KVOPJ(IT),  KVUPJ(IT),  KVDPJ(IT)) 

CALL  UOM(  CC,  AR,  PJM(IM,IL,IG) ,  MOPJM(IT), 

•  MUPJM(IT),  KVOPJM(IT),  KVUPJM(IT),  KVDPJM(IT)) 

CALL  UOM(  CC,  AR,  PJ1(IM,IL,IG) ,  MOPJl(IT), 


MAR07210 
MAR07220 
MAR07230 
MAR07240 
MAR07250 
MAR07260 
MAR07270 
MAR07280 
MAR07290 
MAR07300 
MAR07310 
MAR07320 
MAR07330 
MAR07340 
MAR07350 
MAR07360 
MAR07370 
MAR07380 
MAR07390 
MAR07400 
MAR07410 
MAR07420 
MAR07430 
MAR07440 
MAR07450 
MAR07460 
MAR07470 
MAR07480 
MAR07490 
MAR0750O 
MAR07510 
MAR07520 
MAR07530 
MAR07540 
MAR07550 
MAR07560 
MAR07570 
MAR07580 
MAR07590 
MAR07600 
MAR07610 
MAR07620 
MAR07630 
MAR07640 
MAR07650 
MAR07660 
MAR07670 
MAR07680 
MAR07690 
MAR0770O 
MAR07710 
MAR07720 
MAR07730 
MAR07740 
MAR07750 
MAR07760 


G4 


270 

CONTINUE 

260 

CONTINUE 

250 

CONTINUE 

IF  (KP  .EQ 

FF  =  0 

240 


MUPJl(IT),  KVOPJl(IT),  KVUPJl(IT),  KVDPJl(IT)) 

CALL  UOM(  CC,  AR,  PJ1M( IM, IL, IG) ,  MOPJIM(IT), 

MUPJIM(IT),  KVOPIM(IT),  KVUPIM(IT),  KVBPIM(IT)) 

CALL  UOM(  CC,  AR,  PEB( IM, IL, IG) ,  MOPEB(IT), 

MUPEB(IT),  KVOPEB(IT),  KVUPEB(IT),  KVDPEB(IT)) 

CALL  UOM(  CC,  AR,  PTEB(IM,IL, IG) ,  MOPTEB(IT), 

MUPTEB(IT),  KVOPTB(IT),  KVUPTB(IT),  KVDPTB(IT)) 

CALL  UOM(  CC,  AR,  PEBM( IM,IL, IG) ,  MOPEBM( IT) , 

MUPEBM(IT),  KVOPBM(IT),  KVUPBM( IT) ,  KVDPBM(IT)) 

CALL  UOM(  CC,  AR,  PTEBM(IM,IL,IG) ,  MOPTBM( IT) , 

MUPTBM(IT),  KVOBTM(IT),  KVUBTM(IT),  KVDBTM(IT)) 
ENDIF 


0)  THEN 


ELSE 

FF=FLOAT(KV( IT) )/FLOAT(KP) 
ENDIF 

RRML( IT)=RRML( IT)*FF 
CONTINUE 

COMPUTE  TIME  AVERAGE  MAD 
MEASURE  OF  EFFECTIVENESS 


DO  300  IT=1 


MUPMLE 

MUTSCA 

MUTSAM 

MUPJ 

MUPJM 

MUPJ1 

MUPJ1M 

MUPEB 

MUPTEB 

MUPEBM 

MUPTBM 

MOPMLE 

MOTSCA 

MOTSAM 

MOPJ 

MOPJM 

M0PJ1 

M0PJ1M 

MOPEB 

MOPTEB 

MOPEBM 

MOPTBM 

MDPMLE 

MDTSCA 

MDTSAM 


IT) 

IT) 

IT) 

IT) 

IT) 

IT) 

IT) 

IT) 

(IT) 

IT) 

IT) 

IT) 

IT) 

IT) 

IT) 

IT) 

IT) 

IT) 

IT) 

(IT)  = 

IT) 

IT) 

IT) 

IT) 

IT) 


NVYR 

=  MUPMLE(IT)/KV(IT) 
=  MUTSCA(IT)/KV(IT) 
=  MUTSAM(IT)/KV(IT) 
=  MUPJ  (IT)/KV(IT) 
=  MUPJM(IT)/KV(IT) 
=  MUPJ1(IT)/KV(IT) 
=  MUPJ1M(IT)/KV(IT) 
=  MUPEB(IT)/KV(IT) 

=  MUPTEB(IT)/KV(IT) 
=  MUPEBM(IT)/KV(IT) 
=  MUPTBM(IT)/KV(IT) 
=  MOPMLE(IT)/KV(IT) 
=  MOTSCA(IT)/KV(IT) 
=  MOTSAM(IT)/KV(IT) 
=  MOPJ  (IT)/KV(IT) 
=  MOPJM(IT)/KV(IT) 
=  M0PJ1(IT)/KV(IT) 
=  M0PJ1M(IT)/KV(IT) 
■  MOPEB(IT)/KV(IT) 
MOPTEB(IT)/KV(IT) 
MOPEBM(IT)/KV(IT) 
MOPTBM(IT)/KV(IT) 
MUPMLE(IT)  +  MOPMLE(IT) 
MUTSCA(IT)  +  MOTSCA(IT) 
MUTSAM(IT)  +  MOTSAM(IT) 


MAR07770 
MAR07780 
MAR07790 
MAR07800 
MAR07810 
MAR07820 
MAR07830 
MAR07840 
MAR07850 
MAR07860 
MAR07870 
MAR07880 
MAR07890 
MAR07900 
MAR07910 
MAR07920 
MAR07930 
MAR0  7940 
MAR07950 
MAR07960 
MAR07970 
MAR07980 
MAR07990 
MAR08000 
MAR08010 
MAR08020 
MAR08030 
MAR08040 
MAR08050 
MAR08060 
MAR08070 
MAR08080 
MAR08090 
MAR08100 
MAR08110 
MAR08120 
MAR08130 
MAR08140 
MAR08150 
MAR08160 
MAR08170 
MAR0S180 
MAR08190 
MAR08200 
MAR08210 
MAR08220 
MAR08230 
MAR08240 
MAR08250 
MAR08260 
MAR08270 
MAR08280 
MAR08290 
MAR08300 
MAR08310 
MAR08320 
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MDPJ  (IT)  =  MUPJ  (IT)  + 

MDPJM  (IT)  =  MUPJM  (IT)  + 

MDPJ1  (IT)  =  MUPJ1  (IT)  + 

MDPJIM(IT)  =  MUPJIM(IT)  + 

MDPEB  (IT)  =  MUPEB  (IT)  + 
MDPTEB  (IT)  =  MUPTEB  (IT) 

MDPEBM(IT)  ■  MUPEBM(IT)  + 

MDPTBM(IT)  =  MUPTBM(IT)  + 
300   CONTINUE 

********************************************************************** 

**       PRINT  SECTION  FOR  ALL  OUTPUT.  ** 

***  ft  *  A  ft  ft  ft*********Vafc***************************ft*  ft  *  ft  ****************** 


MOPJ  (IT) 
MOPJM  (IT) 
M0PJ1  (IT) 
MOPJIM(IT) 
MOPEB  (IT) 
+  MOPTEB  (IT) 
MOPEBM(IT) 
MOPTBM( IT) 


** 


905 


** 


906 


907 


PRINT  GROUP  CHARACTERISTICS 
WRITE(6,'(1X,"TEN  OR  LESS  BAYES  VARIANCE  ITERATIONS")') 
WRITE(6,*(1X,I2,"  BAYES  VARIANCE  ITERATION(S)  "  ) ' )  IC 
WRITE(6,'(1X,"XBB  ",F12.6)')  XBB 

T,KK,SSB,SSE,SHRX?SHR,ZB,ZTB,PZTB,A?AT 

/,'  T=',I6,   //  K=',I6, 


/, 


WRITE (6. 9 05) 

FORMAT(f  GROUP  CHARACTERISTICS f 

*  /,'   SSB=' ,F12.4, 

*  /,'   SHR=',F12.4, 

*  /,'  PZTB='  ,F12.4, 
PRINT  RESULTS  FOR 

WRITE(6,*)  '  RESULTS 
WRITE(6,906)  (IVYR(I), 
F0RMAT(14X, 10112) 

WRITE(6.907)  'RTSCA  ',  (RTSCA(IT),  IT=1,NVYR) 
f  ',A14,10F12. 


C  *** 


C 

c 
c 

** 


*** 
*** 
*** 


c   *** 


908 


** 


/,'   SSE=',F12.4, 
/,'    ZB=',F12.4, 

/,'    A=\F12.4,  /,' 
TRANSFORM  SPACE 
FOR  TRANSFORM  SPACE' 
1=1,  NVYR) 


SHRX=' ,F12. 
ZTB=' ,F12.  4, 
AT=" ,F12.  4) 


FORMAT( 
WRITE (6 
WRITE(6 
WRITE  (6 
WRITE(6 
WRITE(6 
WRITE(6 
WRITE(6 
WRITE(6 
WRITE (6 

PRINT 
WRITE(6 
WRITE(6 
WRITE(6 
WRITE(6 
WRITE(6 
WRITE(6 
WRITE(6 
WRITE(6 
WRITE(6 
WRITE(6 
WRITE(6 
WRITE(6 
WRITE(6 
WRITE(6 
WRITE(6 
FORMAT( 
WRITE (6 
WRITE(6 

PRINT 


907) 
907) 
907) 
907) 
907) 
907) 
907) 
907) 
907) 


RESULTS 
i 


*) 

*)  ' 
907) 
907) 
907) 
907) 
907) 
907) 
907) 
907) 
907) 
907) 
907) 
907) 
908) 


RJ 
RJ1 

RXEB 
RXTEB 
RJLT 
RML 
RSL 
RTJ 
RTJLT 
FOR 


4) 


(RJ  (IT), 
(RJ1  (IT), 
(RXEB  (IT), 

(RXTEB  (IT) 
(RJLT  (IT), 
(RML  (IT), 
(RSL  (IT), 
(RTJ  (IT), 
(RTJLT(IT), 
ORIGINAL  SPACE 


IT=1,NVYR) 

IT=1,NVYR) 

IT=1,NVYR) 

,  IT=1,NVYR) 

IT=1,NVYR) 

IT=1,NVYR) 

IT=1,NVYR) 

IT=1,NVYR) 

IT=1,NVYR) 


RESULTS 
RRML 
RRTSCA 
RRTSAM 
RRJ 
RRJM 
RRJ1 
RRJ1M 
RRPEB 
RRPTEB 
RRPEBM 
RRPTBM 
RRJLT 
KV 


FOR  ORIGINAL 
(RRML  (IT), 
(RRTSCA(IT), 
(RRTSAM(IT), 
(RRJ  (IT), 
(RRJM  (IT), 
(RRJ1  (IT), 
(RRJ1M  (IT), 
(RRPEB(IT), 

(RRPTEB(IT) 
(RRPEBM(IT), 
(RRPTBM(IT), 
(RRJLT(IT), 
(KV   (IT), 


,A14, 10112) 


'KP 
t 


SPACE' 
IT=1,NVYR) 

,  IT=1,NVYR) 

IT=1,NVYR) 

IT=1,NVYR) 

IT=1,NVYR) 

IT=1,NVYR) 

IT=1,NVYR) 
IT=1,NVYR) 
,  IT=1,NVYR) 


IT=1,NVYR) 
IT=1,NVYR) 
IT=1,NVYR) 
IT=1,NVYR) 


'(9X, 

*) 

MAD  CALCULATIONS 


,112)')  KP 


MAR08330 
MAR08340 
MAR08350 
MAR08360 
MAR08370 
MAR08380 
MAR08390 
MAR08400 
MAR08410 
MAR08420 
MAR08430 
MAR08440 
MAR08450 
MAR08460 
MAR08470 
MAR08480 
MAR08490 
MAR08500 
MAR03510 
MAR08520 
MAR08530 
MAR08540 
MAR08550 
MAR08560 
MAR08570 
MAR08580 
MAR08590 
MAR08600 
MAR08610 
MAR08620 
MAR08630 
MAR08640 
MAR08650 
MAR08660 
MAR08670 
MAR08680 
MAR08690 
MAR08700 
MAR08710 
MAR08720 
MAR08730 
MAR08740 
MAR08750 
MAR08760 
MAR08770 
MAR08780 
MAR08790 
MAR08800 
MAR08810 
MAR08820 
MAR08830 
MAR08840 
MAR08850 
MAR08860 
MAR08870 
MAR08880 
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WRITE(6,'(2X,' 'PMLE' ')')  MAR08890 
CALL  WRMAD(IVYR,NVYR.MUPMLE,KVUPML,MOPMLE,KVOPML,MDPMLE,KVDPML)   MAR08900 

WRITE(6,'(2X,"TSCA' ')' )  MAR08910 
CALL  WRMAD(IVYR,NVYR  MUTSCA,KVUTSA,MOTSCA,KVOTSArMDTSCAJKVr)TSA)  MAR08920 
WRITE(6,'(2X,* 'TSCAMf ')')  *  MAR08930 
CALL  WRMAD(IVYR,NVYR  MUTSAM,KVUTSM,MOTSAM,KVOTSM,MDTSAM,KVDTSM)   MAR08940 

WRITE(6,'(2X,"PJ")f )  MAR08950 

CALL  WRMAD(IVYR,NVYR,MUPJ,KVUPJ,MOPJ,KVOPJ,MDPJ,KVDPJ)  MAR08960 

WRITE(6,'(2X,"PJM")')  MAR08970 

CALL  WRMAD(IVYR,NVYR,MUPJM,KVUPJM,MOPJM,KVOPJM,MDPJM,KVDPJM)  MAR08980 

WRITE(6,'(2X,"PJl")')  MAR08990 

CALL  WRMAD(IVYR,NVYR.MUPJ1,KVUPJ1,M0PJ1,KV0PJ1JMDPJ1,KVDPJ1)  MAR09000 

WRITE(6,'(2X,"PJ1M' f)')  MAR09010 
CALL  WRMAD(IVYR,NVYRJMUPJ1MJKVUP1M,M0PJ1M,KV0P1M,MDPJ1M,KVDP1M)   MAR09020 

WRITE(6,'(2X,"PEB")')  MAR09030 

CALL  WRMAD(IVYR,NVYR?MUPEB,KVUPEB,MOPEB,KVOPEB,MDPEB,KVDPEB)  MAR09040 

WRITE(6,"(2X,"PEBM' ')')  MAR09050 
CALL  WRMAD(IVYR,NVYR.MUPEBM,KWPBM,MOPEBM,KVOPBM,M)PEBM,KVDPBM)   MAR09060 

WRITE(6,'(2X,"PTEB'  f)')  MAR09070 
CALL  VRMAD(IVYR,NVYR,MUPTEB,KVUPTB,MOPTEB,KVOPTB,MDPTEB,KVDPTB)   MAR09080 

WRITE(6,'(2X,"PTEBM")')  MAR09090 
CALL  WRMAD(IVYR.NVYR,MUPTBM,KWBTM,MOPTBM,KVOBTM,MDPTBM,KVDBTM)   MAR09100 

WRITE(6,*)    '  f  MAR09110 

**       PRINT  THE  CENTRAL  INVENTORY  MATRIX  MAR09120 

VRITE(6,*)  *CINV  MATRIX'  MAR09130 

DO  1006  IT=1.NYR  MAR09140 

WRITE(6,*)  fYEAR',IT  MAR09150 

DO  1005  IL=l,NLOS  MAR09160 

WRITE( 6 , 1007) (CINV( IM, IL, 1 , IT) , IM=1 ,NMOS)  MAR09 170 

1005  CONTINUE  MAR09180 

1006  CONTINUE  MAR09190 

1007  F0RMAT(15(1X,F7.2))  MAR09200 
WRITE(6,*)    '  '  MAR09210 

**       PRINT  THE  ATTRITION  MATRIX  MAR09220 

VRITE(6,*)  'Y  MATRIX1  MAR09230 

DO  1016  IT=1,NYR  MAR09240 

WRITE(6,*)  'YEAR', IT  MAR09250 

DO  1015  IL=l,NLOS  MAR09260 

WRITE( 6 , 1007 ) ( Y( IM, IL, 1 , IT) , IM=1 ,NMOS)  MAR09270 

1015  CONTINUE  MAR09280 

1016  CONTINUE  MAR09290 
WRITE(6,*)    '  '  MAR09300 

**       PRINT  ALL  ARRAYS  MAR09310 

WRITE(6,'(2X,  "ARRAY  TSCA*')')  MAR09320 

CALL  ARRAY(MXMOS,MXLOS,MXGRD,NMOS,NLOS,NGRD,TSCA)  MAR09330 

WRITE(6,'(2X,  "ARRAY  TSCAM")')  MAR09340 

CALL  ARRAY(MXMOS,MXLOS,MXGRD,NMOS,NLOS,NGRD, TSCAM)  MAR09350 

VRITE(6,'(2X,  "ARRAY  SHREB")')  MAR09360 

CALL  ARRAY(MXMOS,MXLOS,MXGRD,NMOS,NLOS,NGRD, SHREB)  MAR09370 

WRITE(6,'(2X,  "ARRAY  SHRTEB")')  MAR09380 

CALL  ARRAY(MXMOS,MXLOS,MXGRD,NMOS,NLOS,NGRD, SHRTEB)  MAR09390 

WRITE(6,'(2X,  "ARRAY  PMLE")')  ,  MAR09400 

CALL  ARRAY(MXMOS,MXLOS,MXGRD,NMOS,NLOS,NGRD,PMLE)  MAR09410 

VRITE(6, '(2X,  "ARRAY  PJ")')  MAR09420 

CALL  ARRAY(MXMOS,MXLOS,MXGRD,NMOS,NLOS,NGRD,PJ)  MAR09430 

WRITE(6,'(2X,  "ARRAY  PJM")1)  MAR09440 
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CALL  ARRAYC MXMOS, MXLOS, MXGRD, NMOS, NLOS, NGRD, PJM) 

WRITE(6,'(2X,  "ARRAY  PJl")') 

CALL  ARRAYC  MXMOS , MXLOS , MXGRD , NMOS , NLOS , NGRD , P Jl ) 

WRITEC  6 , ' ( 2X ,  "  ARRAY  PJ1M "  ) '  ) 

CALL  ARRAYC MXMOS , MXLOS , MXGRD , NMOS , NLOS , NGRD , P JIM) 

WRITE  ( 6,' (2X,"  ARRAY  PEB")') 

CALL  ARRAYC MXMOS , MXLOS , MXGRD , NMOS , NLOS , NGRD , PEB ) 

WRITE(6,'(2X,"ARRAY  PTEB")') 

CALL  ARRAY( MXMOS , MXLOS , MXGRD , NMOS , NLOS , NGRD , PTEB ) 

WRITE(6,'(2X,,'ARRAY  PEBM")') 

CALL  ARRAY( MXMOS, MXLOS, MXGRD, NMOS, NLOS, NGRD, PEBM) 

WRITE(6,,(2X," ARRAY  PTEBM")') 

CALL  ARRAY( MXMOS, MXLOS, MXGRD, NMOS, NLOS, NGRD, PTEBM) 

WRITEC  6,' (2X,"  ARRAY  V"  ) '  ) 

CALL  ARRAY( MXMOS , MXLOS , MXGRD , NMOS , NLOS , NGRD , V) 

WRITE  ( 6,*  (2X,  "ARRAY  XB")') 

CALL  ARRAY( MXMOS , MXLOS , MXGRD , NMOS , NLOS , NGRD , XB ) 

WRITE  ( 6,' C2X,"  ARRAY  VT"  )  *  ) 

CALL  ARRAYC MXMOS, MXLOS, MXGRD, NMOS, NLOS, NGRD, VT) 

WRITEC 6 , ' ( 2X ,  "  ARRAY  XTB  "  ) * ) 

CALL  ARRAY( MXMOS , MXLOS , MXGRD , NMOS , NLOS , NGRD , XTB ) 

WRITE(  6 ,  *  ( 2X,  "  ARRAY  XJ"  ) '  ) 

CALL  ARRAYC MXMOS, MXLOS, MXGRD, NMOS, NLOS, NGRD, XJ) 

WRITE(6,'(2X,' 'ARRAY  XJl")') 

CALL  ARRAY( MXMOS , MXLOS , MXGRD , NMOS , NLOS , NGRD , XJl ) 

WRITE  ( 6,*  (2X,  "ARRAY  XEB")!) 

CALL  ARRAYC  MXMOS , MXLOS , MXGRD , NMOS , NLOS , NGRD , XEB ) 

WRITEC6,'(2X,  "ARRAY  XTEB")') 

CALL  ARRAYC  MXMOS , MXLOS , MXGRD , NMOS , NLOS , NGRD , XTEB ) 
* 

END 

SUBROUTINE  WRMADCIYR,NYR,MU,KU,MO,KO,MD,KD) 

INTEGER  IYRCNYR),KUCNYR),KOCNYR),KDCNYR) 

REAL    MUCNYR),MO(NYR),MDCNYR) 

WRITEC 6,' C/9X,'  'DEVIATION'  '  ,T26,I2,T42,I2,T59,I2/)  '  ) 

*  (IYRCIW),IW=1,NYR) 

WRITEC6,'(9X,'  'UNDERAGE'  '  ,T22,F9.  6,"  C"  ,12,")"  ,T38,F9.6, 

*  "  C",I2,")",T55,F9.6,"  C  '  *  ,12,  "  )  '  '  )  '  ) 

*  MUC1),KUC1),MUC2),KUC2),MU(3),KUC3) 
WRITE(6,'(9X,'  'OVERAGE'  '  ,T22,F9.6,'  '  C'  '  ,12,")"  ,T38,F9.6, 

*  "  C'M2,")",T55,F9.6,"  (", 12, ")")*) 

*  MOC 1) ,KOC 1) ,MOC2) ,KOC2) ,MOC3) ,KOC3) 
WRITE(6,'(9X,"TOTAL",T22,F9.6,"  ("  ,12,")"  ,T38,F9.6, 

*  "  C",I2,")",T55,F9.6,"  C",I2, ")")') 

*  MDC1),KD(1),MDC2),KDC2),MDC3),KDC3) 
WRITEC6,'C/)') 

RETURN 
END 

SUBROUTINE  ARRAYCMM,ML,MG,NM,NL,NG,OBJl) 

REAL  0BJ1CMM,ML,MG) 
DO  10  IG1=1,NG 
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MAR09450 
MAR09460 
MAR09470 
MAR09480 
MAR09490 
MAR09500 
MAR09510 
MAR09520 
MAR09530 
MAR09540 
MAR09550 
MAR09560 
MAR09570 
MAR09580 
MAR09590 
MAR09600 
MAR09610 
MAR09620 
MAR09630 
MAR09640 
MAR09650 
MAR09660 
MAR09670 
MAR09680 
MAR09690 
MAR09700 
MAR09710 
MAR09720 
MAR09730 
MAR09740 
MAR09750 
MAR09760 
MAR09770 
MAR09780 
MAR09790 
MAR09800 
MAR09810 
MAR09820 
MAR09830 
MAR09840 
MAR09850 
MAR09860 
MAR09870 
MAR09880 
MAR09890 
MAR09900 
MAR09910 
MAR09920 
MAR09930 
MAR09940 
MAR09950 
MAR09960 
MAR09970 
MAR09980 
MAR09990 


DO    10    IL1=1.NL 
10  WRITE(6.f(2X,15(F7.3,lX))')    (0BJ1( IM1, IL1,IG1) ,IM1=1,NM) 

WRITE(6,"(/)f) 
RETURN 
END 

SUBROUTINE  FTT(CINV,Y,    X) 

REAL  CINV,X,Y 
TEMP=Y/(1+CINV) 
TEMP1=( 1. +Y)/( 1+CINV) 

TEMP=SQRT(.5+CINV)*. 5*(ASIN( -1. +2.  *TEMP)   + 
*  ASIN(«1.+2.*TEMP1)) 

X=TEMP 
END 

t      t      t     t      i  ._i  ..i__t__t_  .  i . .  i .    i      i      i.    t.    i      i      t-.»    .i-.t.  .*-■!- .t    -t  ..i      t»i      ittttttittttii    .T,,t|    l,,l..f.    Tul    ilut    iliiT.i1..ti    1,  .1  ,.*  I.1..1.  .1.  ,!■  .1.  i  LJt  ill,  ii  ii  ■  !.i.li 

SUBROUTINE  UOM(VINV,ACT,OBJ,OVER,UNDER,KO,KU,KD) 

REAL  OBJ 

IF   (ACT  .  GT.    OBJ)   THEN 

UNDER  =  UNDER+VINV*( ACT-OB J) 

KU=KU+1 
ELSE 

OVER     =  OVER  +VINV*(OBJ-ACT) 

K0=K0+1 
END  IF 
KD=KD+1 
RETURN 
END 

REAL  FUNCTION  RH0(XB,XBB ,SSE,SSB ,DEE,KK,T,NN) 

INTEGER  T 

REAL*8   SIGSQ,    V,    U 

SIGSQ=SSE/(NN+2) 

V=SSB/T 

U=KK*(KK-3)*(XB-XBB)**2/(V*(KK-1)) 

RH0=DMIN1(1.D0,    DEE/DSQRT(U)) 


END 


REAL  FUNCTION  MILLER(    X,    AVGINV   ) 


«■****** 


TM  =  X/SQRT(AVGINV+.  5 )  +  l.  570796 
TEMPI  =  ASIN(SQRT(1. /(AVGINV+1))) 
IF   (TM  .LT.    TEMPI)   THEN 

MILLER  =  0.  0 

GOTO   1 
ENDIF 
IF  (TM  .GT.    3.  14159-TEMP1)  THEN 

MILLER  =1.0 

GOTO   1 
ENDIF 
IF  (COS(TM)    .EQ.    0.0)   THEN 

MILLER  =  0.  5 

GOTO  1 
ENDIF 


MAR10000 
MAR10010 
MAR10020 
MAR10030 
MAR10040 
MAR10050 
MAR10060 
MAR10070 
MAR10080 
MAR10090 
MAR10100 
MAR10110 
MAR10120 
MAR10130 
MARIO  140 
MAR10150 
MAR10160 
MAR10170 
MAR10180 
MAR10190 
MAR10200 
MAR10210 
MAR10220 
MAR10230 
MAR10240 
MAR10250 
MAR10260 
MAR10270 
MAR10280 
MAR10290 
MAR10300 
MAR10310 
MAR10320 
MAR10330 
MAR10340 
MAR10350 
MAR10360 
MAR10370 
MAR10380 
MAR10390 
MAR10400 
MAR10410 
MAR10420 
MAR10430 
MAR10440 
MAR10450 
MAR10460 
MAR10470 
MAR10480 
MAR10490 
MAR10500 
MAR10510 
MAR10520 
MAR10530 
MAR10540 
MAR10550 
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IF   (COS(TM)    .  GT.    0.0)    SIGNUM  =   1. 

IF   (COS(TM)    .LT.    0.0)    SIGNUM  =   -1. 

TEMP2=1-(SIN(TM)+(SIN(TM)-(1.  /SIN(TM)) )/AVGINV)**2 

IF   (TEMP2.LT.  0.0)   TEMP2=0.  0 

MILLER=.  5*( 1.  -SIGNUM*TEMP2**.  5) 

END 


***********•**!■«-*•*■> 


1   ■  I  _>__  I  __l   ■  I    I    I    1    i    t    I    I   .  t  __f  __  I  __  t    t    t    t    I   _!__!_-<    I    I 


REAL  FUNCTION  SCINV(X,A) 

**  SCALE   INVERSION  (TRANSFORM  =>  ORIGINAL  SCALE) 

R=X/SQRT(A+. 5) 
IF(R  .LT.    -1.570796)  THEN 

SCINV=0. 
ELSE   IF(R  .GT.      1.570796)  THEN 

SCINV=1. 
ELSE 

SCINV=.  5*(1.+SIN(R)) 
ENDIF 
END 

REAL  FUNCTION  VAR(    ZZ    ) 

DATA  AA/ 1.6835/,    Bl/-.  8934/,    B2/.  8991/ 
IF  (    ZZ  .GT.    2.2)   THEN 

VAR  =  1. 
ELSE 

IF  (    ZZ 
VAR  = 
ELSE 

VAR  =  AA*(ZZ**B1)*(ZZ-1.  )**B2 
ENDIF 
.    ENDIF 
END 

REAL  FUNCTION  CHISQR( ID,OBJ,INV,ACT) 

REAL  OBJ,INV 

IF   (ID.  NE.  0   .AND.    OBJ.  NE.  0.    .AND.    OBJ. 

*  CHISQR  =  (INV*(ACT-OBJ)**2)/(OBJ*(l. 
END 

SUBROUTINE  READER(MOS, LOS, GRADE, YR,INV,Y,    NMOS.SMOS, 

*  NGRD,SGRD,    SLOSl,SLOS2,    ST1,ST2,    IM,IL,IG,IT,    IEOF) 


LT. 
05 


1.001)   THEN 


NE.  1.  ) 
-OBJ)) 


READ 
INTEGER 
REAL  Y 
INTEGER 
INTEGER 


RECORD  AND  ACCUMULATE   LOSSES 

MOS,    LOS,    GRADE,    YR,    INV,      LOSS(8) 


1 
100 


ST1,    ST2 

SLOSl,    SL0S2 
INTEGER  SM0S(40),    NMOS 
INTEGER  SGRD(IO),    NGRD 
ICNT=0 

READ( 1,100, END=999 )   YR , MOS , GRADE , LOS , INV , ( L0SS( I ) , 1=1 , 8) 
FORMAT(4I2,9I5) 

CHECK  IF  RECORD  MEETS  SELECTION  CRITERIA.      OTHERWISE  REJECT. 
IM=0 


MAR10560 
MAR 105 70 
MAR10580 
MAR10590 
MAR10600 
MAR10610 
MAR10620 
MAP.10630 
MAR10640 
MAR10650 
MAR10660 
MAR10670 
MAR10680 
MAR10690 
MAR10700 
MAR10710 
MAR10720 
MAR10730 
MAR10740 
MAR10750 
MAR10760 
MAR10770 
MAR10780 
MAR10790 
MAR10800 
MAR10810 
MAR10820 
MAR10830 
MAR 10 840 
MAR10850 
MAR10860 
MAR10870 
MAR10880 
MAR10890 
MAR10900 
MAR10910 
MAR10920 
MAR10930 
MAR10940 
MAR10950 
MAR10960 
MAR10970 
MAR10980 
MAR10990 
MAR11000 
MAR11010 
MAR11020 
MAR11030 
MAR11040 
MAR11050 
MAR11060 
MAR11070 
MAR11080 
MAR11090 
MAR11100 
MAR11110 
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10 

* 
20 


30 
40 


tHt 


** 


50 
999 


DO  10  I=1,NM0S 
IF(M0S  .EQ.  SMOS(I))  THEN 
IM=I 

GO  TO  20 
END  IF 
CONTINUE 
GO  TO  1 

CONTINUE 
DO  30  I=1,NGRD 
IF( GRADE   .EQ.    SGRD(I))  THEN 
IG=I 

GO  TO  40 
END  IF 
CONTINUE 
GO  TO   1 
CONTINUE 

IF(L0S   .  LT.    SL0S1   .OR.    LOS   .  GT.    SL0S2)   GO  TO   1 
IL--L0S-SL0S1+1 

YEARS  OVER  ST2  ARE  USED  FOR  VALIDATION 
IF(YR  .LT.    ST1)   GO  TO   1 
IT=YR-ST1+1 

COMPUTE  TOTAL  LOSS 
Y=0 
DO  50  1=1,8 

Y=Y  +  LOSS(I) 
CONTINUE 

RETURN 
IE0F=1 
END 


MAR1 
MAR1 
MAR1 
MAR1 
MAR1 
MAR1 
MAR1 
MAR1 
MAR1 
MAR1 
MAR1 
MAR1 
MAR1 
MAR1 
MAR1 
MAR1 
MAR1 
MAR1 
MAR1 
MAR1 
MAR1 
MAR1 
MAR1 
MAR1 
MAR1 
MAR1 
MAR1 
NAR1 
MAR1 
MAR1 
MAR1 
MAR1 


1120 
1130 
1140 
1150 
1160 
1170 
1180 
1190 
1200 
1210 
1220 
1230 
1240 
1250 
1260 
1270 
1280 
1290 
1300 
1310 
1320 
1330 
1340 
1350 
1360 
1370 
1380 
1390 
1400 
1410 
1420 
1430 
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APPENDIX  C.     SIMULATION  FOR  CHOICE  OF  AVERAGE 


INVENTORY  VALUES 


PROGRAM  MCSIM 


SIMULATION  COMPARISON  OF  INVENTORY  MEAN  PERFORMANCE 

PURPOSE:    TO  IDENTIFY  THE  BEST  MEAN  VALUE  (ARITHMETIC, 

GEOMETRIC  AND  HARMONIC)   FOR  USE   IN  THE   INVERSION  OF 
TRANSFORMED  DATA  TO  THE  ORIGINAL  SPACE.   ' 

PROGRAMMER:  CAPT  C  R  DICKINSON,  USMC 

VARIABLE  DESCRIPTION: 

LAMBDA  -  PARAMETER  USED  IN  THE  POISSON  RANDOM  NUMBER 

GENERATOR 
TRIALS  -  VECTOR  OF  INVENTORY  VALUES  GENERATED  FROM 

THE  POISSON  RANDOM  NUMBER  GENERATOR 
REPS   -  NUMBER  OF  INVENTORY  VALUES  TO  BE  GENERATED 
BASIC   -  BASIC  INVERSION  OF  TRANSFORMED  VALUES 
FTE1        -  FREEMAN-TUKEY  EXACT  INVERSION  WITH 


FTE2 


FREEMAN-TUKEY  EXACT  INVERSION  WITH 


******* 

■k 
•k 
■k 

* 

* 

* 

■k 
■k 
* 
* 
■k 
■k 
■k 
■k 
■k 
■k 
k 
•k 
■k 


Z  -  TRANSFORMED  SCALE  VARIABLE 

MEANS   -  VECTOR  OF  ARITHMETIC,  GEOMETRIC  AND  HARMONIC 

MEANS 
TRSUM  -  SUMMATION  OF  TRIALS  FOR  USE  IN  THE 

ARITHMETIC  MEAN 
TRMULT  -  MULTIPLICATION  OF  TRIALS  FOR  USE  IN  THE 

GEOMETRIC  MEAN 
TRINV  -  SUMMATION  OF  THE  INVERSE  OF  TRIALS  FOR  USE 

IN  THE  HARMONIC  MEAN 

INPUT/OUTPUT:  * 

INPUT  -  HARD  CODED  PARAMETERS  * 

OUTPUT  -  SIM  OUT  Al  (DISK  FILE  03)  * 

**-k-k-kirk-k*-kirk-k-k-!rt<*-irt:-k-k-k-trk-k-k-k-k-k  A  h  i'c  k  k******irk*1rkirk*itJrkirirteir*it******ir*-k-k-k-k-k 
•k 

***      VARIABLE  DECLARATION 

INTEGER  TRIALS( 100) ,Y( 100) ,REPS(3) ,PROB,R, 
C     COUNT, TRSUM, YTEMP,YSUM 

REAL  Z( 100) ,T( 100) ,MEAN(3) ,BASIC(3) ,DIFF(3) ,LAMBDA(5) , 
C     P(4),FTE1(3),FTE2(3),TRMULT 

DOUBLE  PRECISION  DPSEED 


* 

VtAtV 


CALL  EXCMS('FILEDEF  02  DISK  REPLIC  LISTING  Al') 

VARIABLE  INITIALIZATION 
PI=3.  14159 
DPSEED=889246.D0 
LAM=5 
R=3 

PR0B=4 
COUNT=2 


MCS00010 
MCS00020 
MCS00030 
MCS00040 
MCS00050 
MCS00060 
MCS00070 
MCS00080 
MCS00090 
MCS00100 
MCS00110 
MCS00120 
MCS00130 
MCS00140 
MCS00150 
MCS00160 
MCS00170 
MCS00180 
MCS00190 
MCS00200 
MCS00210 
MCS00220 
MCS00230 
MCS00240 
MCS00250 
MCS00260 
MCS00270 
MCS00280 
MCS00290 
MCS00300 
MCS00310 
MCS00320 
MCS00330 
MCS00340 
MCS00350 
MCS00360 
MCS00370 
MCS00380 
MCS00390 
MCS00400 
MCS00410 
MCS00420 
MCS00430 
MCS00440 
MCS00450 
MCS00460 
MCS00470 
MCS00480 
MCS00490 
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DATA  REPS/20,40,60/ 

DATA  LAMBDA/ 1.  ,2.  ,4.  ,8.  ,16./ 

DATA  P/.05,.  1,.2,.4/ 


* 


10 


tHV* 


HEADER  PRINT 
WRITE  (02,'(/T6,' 'LAMBDA' ' ,T15,' 'REPS' ' ,T27,' 'P"  , 
C      T34,' 'METHOD' *  ,T46  ' 'ARITHMETIC' ' ,T59, 
C      "GEOMETRIC"  ,T71,   HARMONIC"/)') 

LOOP  FOR  LAMBDA  PARAMETER 
DO  500  1=1, LAM 
DO  450  J=1,R 

CALL  GGPOS(LAMBDA(I),DPSEED,REPS(J),TRIALS,IER) 

NSUM=0 

DO  10  J1=1,REPS(J) 

TRIALS(J1)=TRIALS(J1)+1 

NSUM=NSUM+TRIALS(J1) 


LOOP  FOR  PROBABILITY  PARAMETER 
DO  400  K=l,PROB 
C0UNT=C0UNT+4 
LINE  COUNTER  FOR  PAGE  BREAK 
IF  (COUNT. GT. 55)  THEN 
WRITE  (2,'(1H1)') 
COUNT=2 

WRITE  ( 02,' (/T6,"  LAMBDA" 
T34,' 'METHOD' ' 


C 
C 


•T15,' 'REPS' ' ,T27,' 'P*  * 
,T46    ' 'ARITHMETIC"  ,T59, 
"GEOMETRIC"  ,T71,      HARMONIC'  '/)') 


c 
c 


20 


it     i  _ 
II  p*  w€ 


ENDIF 

TRSUM=0 

TRINV=0 

TRMULT=0. 

ARCSUM=0. 

TSUM=0. 

YSUM=0 

LOOP  FOR  REPITITION  PARAMETER 
DO  20  K1=1,REPS(J) 

CALL  GGBN(DPSEED,1, TRIALS (K1),P(K),YTEMP) 

YSUM=YSUM+YTEMP 

TRSUM=TRSUM+TRIALS(K1) 

TRMULT=TRMULT+ALOG( REAL( TRIALS ( Kl )) ) 

TRINV=TRINV+1.  /REAL(TRIALS(K1)) 

Y(K1)=YTEMP 

Z(K1)=.  5*(TRIALS(K1)+.  5)**.  5*(ASIN(2.  *Y(K1)/ 

(TRIALS(K1)+1.  )-l.  )+ASIN(2.*(Y(Kl)+l.  )/ 

(TRIALS(K1)  +  1.  )-l.  )+PI) 
T(K1)=Z(K1)/((TRIALS(K1)+.  5)**. 5) 
ARCSUM=ARCSUM+Z(K1) 
TSUM=TSUM+T(K1) 

ZAVG=ARCSUM/REAL(REPS( J) ) 
TAVG=TSUM/REAL(REPS( J) ) 
MEAN  VALUE  CALCULATION 

••''•■* 
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MCS00500 
MCS00510 
MCS00520 
MCS005  30 
'MCS00540 
MCS00550 
MCS00560 
MCS00570 
MCS00580 
MCS00590 
MCS00600 
MCS00610 
MCS00620 
MCS00630 
MCS00640 
MCS00650 
MCS00660 
MCS00670 
MCS00680 
MCS00690 
MCS00700 
MCS00710 
MCS00720 
MCS00730 
MCS00740 
MCS00750 
MCS00760 
MCS00770 
MCS00780 
MCS00790 
MCS00800 
MCS00810 
MCS00820 
MCS00830 
MCS00840 
MCS00850 
MCS00860 
MCS00870 
MCS00880 
MCS00890 
MCS00900 
MCS00910 
MCS00920 
MCS00930 
MCS00940 
MCS00950 
MCS00960 
MCS00970 
MCS00980 
MCS00990 
MCS01000 
MCS01010 
MCS01020 
MCS01030 
MCS01040 
MCS01050 


MEAN(1)=TRSUM/REAL(REPS(J)) 
MEAN(2)=EXP(TRMULT*(1./REAL(REPS(J)))) 
MEAN( 3)=REAL(REPS( J) )/TRINV 

IF  (COS(TAVG).GE.  0.0)  THEN 

SIGNUM=1. 
ELSE 

SIGNUM=-1. 
ENDIF 

DO  30  K2=l  3 

TAVG1=ZAVG/(MEAN(K2)+.  5)**. 5 
IF  (COS(TAVGl).GE.  0.0)  THEN 

SIG1=1. 
ELSE 

SIG1=-1. 
ENDIF 
***      BASIC  AND  FTE  INVERSE  METHOD  CALCULATION 

FTE1(K2)=.  5*(1.  -(SIGNUM*(1.  -( ABS(SIN(TAVG)  + 
C  (SIN(TAVG)-(l.  /SIN(TAVG)))/MEAN(K2)))**2.  )**.  5)) 

FTE2(K2)=.  5*(1.  -(SIG1*(1.  -(ABS(SIN(TAVG1)  + 
C  (SIN(TAVGl)-(l.  /SIN(TAVG1)))/MEAN(K2)))**2.  )**.  5)) 

BASIC(K2)=.  5*(1.+SIN(ZAVG/((MEAN(K2)+.  5)**.  5) 
C  -(PI/2.))) 

OUTPUT  STATEMENTS 
DO  40  K3=l,3 

IF  (K3.EQ.  1)  THEN 

WRITE  (02,'(T6,F4.  1,T15  ,I3,T25  ,F5.  3,T35,  "  BASIC"  , 
C  T46,F7.4,T59,F7.4,T71,F7.4)') 

C  LAMBDA(I),REPS(J),P(K),(BASIC(IK),IK=1,3) 

ELSE 

IF  (K3.EQ.  2)  THEN 

WRITE  (02, '(T35,1  *FTEl' ' ,T46,F7. 4,T59 ,F7.  4, 
C  T71,F7.4)f)  (FTE1(IK),IK=1,3) 

ELSE 

WRITE  (02,*(T35,"FTE2f  *  ,T46,F7.  4,T59,F7.  4, 
C  T71.F7.4/)')  (FTE2(IK),IK=1,3) 

ENDIF 
ENDIF 


* 

30 


40 

CONTINUE 

* 

400 

CONTINUE 

450 

CONTINUE 

500 

CCNTINUE 

STOP 

END 

MCS01060 
MCS01070 
MCS01080 
MCS01090 
MCS01100 
MCS01110 
MCS01120 
MCS01130 
MCS01140 
MCS01150 
MCS01160 
MCS01170 
MCS01180 
MCS01I90 
MCS01200 
MCS01210 
MCS01220 
MCS01230 
MCS01240 
MCS01250 
MCS01260 
MCS01270 
MCS01280 
MCS01290 
MCS01300 
MCS01310 
MCS01320 
MCS01330 
MCS01340 
MCS01350 
MCS01360 
MCS01370 
MCS01380 
MCS01390 
MCS01400 
MCS01410 
MCS01420 
MCS01430 
MCS01440 
MCS01450 
MCS01460 
MCS01470 
MCS01480 
MCS01490 
MCS01500 
MCS01510 
MCS01520 
MCS01530 
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